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PRINCIPAL NOTATION 


SYMBOLS 


Unless specified otherwise, all variables are nondimensional. 
Symbol Definition 


A, B, C 

A', B', C' 

Cpj C v 

E, F 

E, F 

Er 

Ek, Fk 
Ek, Fk 

Eki, Fkj 
Ey v F y 2 

F, G 

A A 

F, G 

hr 

H, H k 

A A 

h,h k 

ij 

J 

k 

k 

k h k t 

Lr 

N eq 


Speed of sound. 

Coefficient submatrices in block tridiagonal system of equations. 

Coefficient submatrices for boundary conditions. 

Specific heats at constant pressure and volume. 

Inviscid flux vectors in the Cartesian or cylindrical coordinate form of the govern- 
ing equations. 

Inviscid flux vectors in the computational coordinate form of the governing 
equations. 

Total energy per unit volume. 

Viscous flux vectors in the Cartesian or cylindrical coordinate form of the govern- 
ing equations. 

Viscous flux vectors in the computational coordinate form of the governing 
equations. 

Non-cross derivative viscous flux vectors in the computational coordinate form of 
the governing equations. 

Cross derivative viscous flux vectors in the computational coordinate form of the 
governing equations. 

Flux vectors in the Cartesian or cylindrical coordinate form of the k-z turbulence 
model equations. 

Flux vectors in the computational coordinate foim of the k-z turbulence model 
equations. 

Stagnation enthalpy per unit mass. 

Non-derivative inviscid and viscous terms in the Cartesian coordinate form of the 
governing equations for axisymmetric flow. 

Non-derivative inviscid and viscous terms in the computational coordinate form 
of the governing equations for axisymmetric flow. 

Grid indices in the £ and rj directions. 

Jacobian matrix of the generalized grid transformation. 

Effective thermal conductivity coefficient. 

Turbulent kinetic energy. 

Laminar and turbulent thermal conductivity coefficient. 

Dimensional reference length. 

Number of governing equations being solved. 
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Svmbol 


Definition 


N u N 2 

P 

Pr r 

Pn, Pr, 


Qx, <]r 

<Jx, q y 

Q 


A 

Q 

R 

Re r 

S 

S' 

S, T 


S, T 

t 

T 

u, v 
u, v, w 

W 


A 

w 

x, r 


x,y 


Number of grid points in the £ and rj directions. 

Static pressure. 

Reference Prandtl number. 

Laminar and turbulent Prandtl number. 

Heat fluxes in the cylindrical x and r directions. 

Heat fluxes in the Cartesian x and y directions. 

Vector of dependent variables in the Cartesian or cylindrical coordinate form of the 
governing equations. 

Vector of dependent variables in the computational coordinate form of the gov- 
erning equations. 

Gas constant. 

Reference Reynolds number. 

Source term subvector in block tridiagonal system of equations. 

Source term subvector for boundary conditions. 

Non-derivative terms in the Cartesian or cylindrical coordinate form of the k-z 
turbulence model equations. 

Non-derivative terms in the computational coordinate form of the k-z turbulence 
model equations. 

Physical time. 

Static temperature. 

Velocities in the Cartesian x and y directions. 

Velocities in the cylindrical x, r, and swirl directions. 

Vector of dependent variables in the Cartesian or cylindrical coordinate form of the 
k-z turbulence model equations. 

Vector of dependent variables in the computational coordinate form of the k-z 
turbulence model equations. 

Cylindrical axial and radial coordinates. 

Cartesian coordinates. 


y 

5 

A, V 

z 

tf 

£/ 

zf\ etc. 
# 1 , 02 , 03 

*2, K 4 

JL 


Ratio of specific heats, c P lc v . 

Difference operator. 

First-order forward and backward difference operators. 

Turbulent dissipation rate. 

Second- and fourth-order explicit artificial viscosity coefficients in constant coeffi- 
cient model. 

Implicit artificial viscosity coefficient. 

Second- and fourth-order artificial viscosity coefficients in nonlinear coefficient 
model. 

Parameters determining type of time differencing used. 

Constants in nonlinear coefficient artificial viscosity model. 

Effective second coefficient of viscosity. 
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Symbol 


Definition 


h, h 

P 

Ph Pt 




p 

o 


Txxi yi etc. 


Laminar and turbulent second coefficient of viscosity. 

Effective viscosity coefficient. 

Laminar and turbulent viscosity coefficient. 

Laminar kinematic viscosity. 

Computational coordinate directions. 

Static density. 

Pressure gradient scaling parameter in nonlinear coefficient artificial viscosity 
model. 

Computational time. 

Elements of shear stress tensor. 

Spectral radius in nonlinear coefficient artificial viscosity model. 


SUBSCRIPTS 

Subscript Definition 


/, j Denotes grid location in £ and r\ directions. 

r Denotes dimensional reference condition. 

t Denotes differentiation with respect to physical time. 

x, r Denotes differentiation with respect to cylindrical coordinate directions. 

jc , y Denotes differentiation with respect to Cartesian coordinate directions. 

rj Denotes differentiation with respect to computational coordinate directions. 

t Denotes differentiation with respect to computational time. 


SUPERSCRIPTS 

Superscript Definition 

n Denotes time level. 

* Denotes solution after first ADI sweep. 
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PROTEUS TWO-DIMENSIONAL 
NAVIER-STOKES COMPUTER CODE - VERSION 2.0 

Volume 1 - Analysis Description 

Charles E. Towne, John R. Schwab, Trong T. Bui 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 


SUMMARY 

A computer code called Proteus has been developed to solve the two-dimensional planar or 
axisymmetric, Reynolds-averaged, unsteady compressible Navier-Stokes equations in strong conservation 
law form. The objective in this effort has been to develop a code for aerospace propulsion applications that 
is easy to use and easy to modify'. Code readability, modularity, and documentation have been emphasized. 

The governing equations are written in Cartesian coordinates and transformed into generalized 
nonorthogonal body-fitted coordinates. They are solved by marching in time using a fully-coupled 
alternating-direction-implicit solution procedure with generalized first- or second-order time differencing. 
The boundary conditions are also treated implicitly, -and may be steady or unsteady. Spatially periodic 
boundary conditions are also available. All terms, including the difhision terms, are linearized using 
second-order Taylor series expansions. Turbulence is modeled using either an algebraic or two-equation 
eddy viscosity model. 

The program contains many operating options. The governing equations may be solved for two- 
dimensional planar flow, or axisymmetric flow with or without swirl. The thin-layer or Euler equations 
may be solved as subsets of the Navier-Stokes equations. The energy equation may be eliminated by the 
assumption of constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- 
and post-shock oscillations in supersonic flow and to minimize odd-even decoupling caused by central 
spatial differencing of the convective terms in high Reynolds number flow. Several time step options are 
available for convergence acceleration, including a locally variable time step and global time step cycling. 
Simple Cartesian or polar grids may be generated internally by the program. More complex geometries 
require an externally generated computational coordinate system. 

The documentation is divided into three volumes. Volume 1, the current volume, is the Analysis De- 
scription, and presents the equations and solution procedure used in Proteus. It describes in detail the 
governing equations, the turbulence model, the linearization of the equations and boundary conditions, the 
time and space differencing formulas, the ADI solution procedure, and the artificial viscosity models. Vol- 
ume 2 is the User's Guide, and contains information needed to run the program. It describes the program's 
general features, the input and output, the procedure for setting up initial conditions, the computer resource 
requirements, the diagnostic messages that may be generated, the job control language used to run the 
program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed informa- 
tion useful when modifying the program. It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram. 

Version 1.0 of the two-dimensional Proteus code was released in late 1989. The current documentation 
covers Version 2.0, released in early 1992. 
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1.0 INTRODUCTION 


Much of the effort in applied computational fluid dynamics consists of modifying an existing program 
for whatever geometries and flow regimes are of current interest to the researcher. Unfortunately, nearly 
all of the available non-proprietary programs were started as research projects with the emphasis on dem- 
onstrating the numerical algorithm rather than ease of use or ease of modification. The developers usually 
intend to clean up and formally document the program, but the immediate need to extend it to new ge- 
ometries and flow regimes takes precedence. 

The result is often a haphazard collection of poorly written code without any consistent structure. An 
extensively modified program may not even perform as expected under certain combinations of operating 
options. Each new user must invest considerable time and effort in attempting to understand the underlying 
structure of the program if intending to do anything more than run standard test cases with it. The user s 
subsequent modifications further obscure the program structure and therefore make it even more difficult 
for others to understand. 

The Proteus two-dimensional Navier-Stokes computer program is a user-oriented and easily-modifiable 
flow analysis program for aerospace propulsion applications. Readability, modularity, and documentation 
were primary objectives during its development. The entire program was specified, designed, and imple- 
mented in a controlled, systematic manner. Strict programming standards were enforced by immediate peer 
review of code modules; Kemighan and Plauger (1978) provided many useful ideas about consistent pro- 
gramming style. Every subroutine contains an extensive comment section describing the purpose, input 
variables, output variables, and calling sequence of the subroutine. With just three clearly-defined ex- 
ceptions, the entire program is written in ANSI standard Fortran 77 to enhance portability. A master ver- 
sion of the program is maintained and periodically updated with corrections, as well as extensions of general 
interest (e.g., turbulence models.) 

The Proteus program solves the unsteady, compressible, Reynolds-averaged Navier-Stokes equations in 
strong conservation law form. The governing equations are written in Cartesian coordinates and trans- 
formed into generalized nonorthogonal body-fitted coordinates. They are solved by marching in time using 
a fully-coupled altemating-direction-implicit (ADI) scheme with generalized time and space differencing 
(Briley and McDonald, 1977; Beam and Warming, 1978). Turbulence is modeled using either the Baldwin 
and Lomax (1978) algebraic eddy- viscosity model or the Chien (1982) two-equation model. All terms, in- 
cluding the diffusion terms, are linearized using second-order Taylor series expansions. The boundary 
conditions are treated implicitly, and may be steady or unsteady. Spatially periodic boundary conditions 
are also available. 

The program contains many operating options. The governing equations may be solved for two- 
dimensional planar flow, or axisymmetric flow with or without swirl. The thin-layer or Euler equations 
may be solved as subsets of the Navier-Stokes equations. The energy equation may be eliminated by the 
assumption of constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- 
and post-shock oscillations in supersonic flow and to minimize odd-even decoupling caused by central 
spatial differencing of the convective terms in high Reynolds number flow. Several time step options are 
available for convergence acceleration, including a locally variable time step and globed time step cycling. 
Simple grids may be generated internally by the program; more complex geometries require external grid 
generation, such as that developed by Chen and Schwab (1988). 

The documentation is divided into three volumes. Volume 1, the current volume, is the Analysis De- 
scription, and presents the equations and solution procedure used in Proteus. It describes in detail the 
governing equations, the turbulence model, the linearization of the equations and boundary conditions, the 
time and space differencing formulas, the ADI solution procedure, and the artificial viscosity models. Vol- 
ume 2 is the User's Guide, and contains information needed to run the program. It describes the program's 
general features, the input and output, the procedure for setting up initial conditions, the computer resource 
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requirements, the diagnostic messages that may be generated, the job control language used to run the 
program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed informa- 
tion useful when modifying the program. It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram. 

Version 1.0 of the two-dimensional Proteus code was released in late 1989 (Towne, Schwab, Benson, 
and Suresh, 1990). The current documentation covers Version 2.0, released in early 1992. 

The authors would like to acknowledge the significant contributions made by their co-workers. Tom 
Benson provided part of the original impetus for the development of Proteus, and did the original coding 
of the block tri-diagonal inversion routines. Simon Chen did the original coding of the Baldwin- Lomax 
turbulence model, and consulted in the implementation of the nonlinear coefficient artificial viscosity model. 
W illiam K nnik developed the original code for computing the metrics of the generalized nonorthogonal grid 
transformation. Fr ank Molls has created separate diagonalized and patched-grid versions of the code. 
Ambady Suresh did the original coding for the second-order time differencing and for the nonlinear coeffi- 
cient artificial viscosity model. These people, along with Dick Cavicchi, Julie Conley, Jason Solbeck, and 
Pat Zeman, have also run many debugging and verification cases. 
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2.0 GOVERNING EQUATIONS 


2.1 GOVERNING EQUATIONS IN CARTESIAN COORDINATES 

The basic governing equations are the two-dimensional compressible Navier-Stokes equations. These 
equations may be found in several standard references (e.g., Hushes and Gaylord, 1964; Schlichting, 1968; 
White, 1974; Anderson, Tannehill, and Pletcher, 1984). In Cartesian coordinates, the two-dimensional 
planar equations 1 can be written in strong conservation law form using vector notation as 

5 Q . dL-HJbL . ^Ll oh 

dt dx dy dx dy 


where 


Q = [p pw pv E r ] r 


(2.2a) 


E = 


pu 

pi? +p 

puv 

( E t + p)u 


pv 

puv 

2 . 

pv +p 
(E r + P)v 


(2.2b) 


(2.2c) 


E„ = 


Re r 


(2.2d) 


** 1 

m xx + VT xy - -frT 4x 




xy 

Tyy 


(2.2e) 


^xy ^yy 


Pr r q y 


Equation (2.1) thus represents, in order, the continuity, x-momentum, ^-momentum, and energy equations, 
with dependent variables p, pu, pv, and E r . 


1 Proteus can be used for both two-dimensional planar or axisymmetric flow. However, the axisymmetric equations 
have some additional terms that complicate the analysis somewhat. For the sake of clarity, the main body of this 
report describes the two-dimensional planar analysis, and the axisymmetric analysis is described in Appendix B. 
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The shear stresses and heat fluxes are given by 



In these equations, t represents time; x and y represent the Cartesian coordinate directions; u and v are 
the velocities in the x and y directions; p, p, and T are the static density, pressure, and temperature; Et is 
the total energy per unit volume; and p, A, and k are the coefficient of viscosity, second coefficient of 
viscosity, and coefficient of thermal conductivity. 

All of the above equations have been nondimensionalized using appropriate normalizing conditions. 
Lengths have been nondimensionalized by L, , velocities by Ur, density by p n temperature by T n viscosity 
by fi n thermal conductivity by k r , pressure and total energy by p r u}, and time by Lr\u r . The reference 
Reynolds and Prandtl numbers are thus defined as Re r = p r u r Lrjp r and Pr, = p r i$jk r T r . 2 * 

Turbulence is modeled using the Boussinesq approach (Schlichting, 1968). The equations presented in 
this section are thus used for both laminar and turbulent flow. For turbulent flow they represent the 
Reynolds time-averaged form of the Navier-Stokes equations, with density fluctuations neglected. They 
may also be interpreted as the Favre or mass-weighted time-averaged form of the equations. With Favre 
time averaging, however, the velocities and thermal variables represent mass-averaged quantities defined by 
u = ~puj~p, etc., where the overbar represents a conventional Reynolds time-averaged quantity. Details on 
Reynolds and Favre time-averaging procedures may be found in Cebeci and Smith (1974), and in Anderson, 
Tannehill, and Pletcher (1984). In either case, p, A, and k represent effective coefficients. For example, in 
turbulent flow pi = pi + p r , where p t and p t are the laminar and turbulent viscosity coefficients, and p t comes 
from some appropriate turbulence model. The models currently available in the Proteus code are the al- 
gebraic eddy viscosity model of Baldwin and Lomax (1978) and the two-equation model of Chien (1982), 
implemented as described in Section 9.0. 

2.2 EQUATION OF STATE 

In addition to the equations presented above, an equation of state is required to relate pressure to the 
dependent variables. Any appropriate equation, or even table, could be used. The equation currently built 
into the Proteus code is the equation of state for thermally perfect gases, p= pRT, where R is the gas con- 
stant. For calorically perfect gases, this can be rewritten as 

p=(y- p ( u1 + y2 )] ( 2 - 4 ) 

where y is the ratio of specific heats, q,/c*. Here the gas constant and specific heats have been 
nondimensionalized by u?jT r . 

If the flow is such that we can assume a perfect gas with constant stagnation enthalpy, the energy 
equation may be eliminated. This assumption is reasonable, for example, in inviscid regions, and in 


2 Note that this Prandtl number does not have a physically meaningful value, but is merely defined by a combination 

of the normalizing conditions for c p , p, and k that appear when the equations are nondimensionalized. 
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adiabatic wall boundary' layers if the Prandtl number is near 1 (Briley and McDonald, 1977). The stag- 
nation enthalpy is defined as 


h T =c p T+-j ( u 2 + v 2 ) (2.5) 

Here the stagnation enthalpy is nondimensionalized by u}. The temperature is thus 

7-=^[v-y ( “ 2 + v 2) ] (2.6) 

and the equation of state becomes 

p = - y 1 p[ A r-y (“ 2 + y2 )] ( 2 - 7 ) 

This equation of state does not require the total energy £r, and the energy equation need not be solved. 
The total energy may be computed from 

£y’ == p/zy — p (2*^) 


2.3 GENERALIZED GRID TRANSFORMATION 


Because the governing equations in the previous section are written in Cartesian coordinates, they are 
not well suited for general geometric configurations. For most applications a body-fitted coordinate system 
is desired. This greatly simplifies the application of boundary conditions and the bookkeeping in the nu- 
merical method used to solve the equations. The following generalized grid transformation, which can be 
orthogonal or nonorthogonal, is therefore used to transform the governing equations from physical (x 7 y, t) 
coordinates to rectangular orthogonal computational (£, rj 7 t ) coordinates. 

£ = £(x 7 y 7 t) 

r, = >/(*, .M) (2-9) 

i — t 

In Proteus, the spatial computational domain is square, with £ and rj each running from 0 to 1. Using the 
chain rule for partial differentiation, the derivatives in the Cartesian foim of the governing equations can 
be replaced using the following expressions. 


_d_ 

dx 

_d_ 

dy 

d_ 

dt 




JL 

d£ + Vx dr. 


-P J- + „ A. 


dt 

d 


dr, 

d 


V v V/ , V 

~ ~~ d£ + + 


drj 


d_ 

dr 


(2.10) 


In the above equations, and in those to follow, subscripts jc and y, or £ and rj 7 denote partial differentiation 
in that coordinate direction. The only task remaining, then, is to develop expressions for the metric coef- 
ficients £ x , rj X7 etc. In differential form we can write 


d£ = £ x dx + Z y tfy + Zjdt 
drj = r\ x dx + *l y dy 4- r\jdt 
dr = dt 


In matrix form this becomes 
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'<«] \i x ty • £,!<&' 

dr] = r) x r\ y tj, dy 

di J [O 0 lj[d/ 

Similarly, 

dbc x T " 

<=fy = J>,, J> T drj 

dt\ Lo 0 1 jl* 

Therefore, 

'ix iy £rl f*C x v x z 
Vx Vy Vt ~ J>J A J>t 

0 0 1J LO 0 1. 

After taking the inverse, 


'i X 

«/ 

it' 

\ y* 


1 

J 

1 

J 

Vx 

Vy 

Vt 

= J -y { 

x s 

y^ x z~ v* 

0 

0 

1 

L 0 

0 

i \J . 


where J is the Jacobian of the transformation, 

J= d(g.n) _ ix Zy 

d(xj>) Vx r, y 

J = Z x Vy - Z y v x 

This can be evaluated from the known physical (x. y) coordinates by noting J = 1 //” 1 and 

! = d(xjQ = x ? x,, 

y$ yr, 

J~ 1 = X ^r, ~ V5 

The metric coefficients themselves are 

£* = J y v 

Zy = ~ Jx n 

v x = ~ J yi. 

Vy = JX{ 

z t =- x Jx-y^y 

V,= ~ x,r, x - y^ y 


( 2 . 11 ) 


( 2 . 12 ) 


( 2 . 13 ) 
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Unless the physical coordinates (x, y) are defined analytically as functions of the computational coordinates 
(<f, rj), the metric coefficients must be computed numerically. 

2.4 GOVERNING EQUATIONS IN COMPUTATIONAL COORDINATES 


Applying the generalized grid transformation of the previous section to equation (2.1) yields 

Q t + + Q v *h + E^ x 4- E T] r\ x + F fty + F v rj y - E v £ x - E v r\ x - F v £ y - F v r\ y = 0 (2. 14) 


This equation is in chain-rule, or weakly conservative form. That is, the conservation flow variables are 
used, but the metrics appear as coefficients of the derivatives instead of inside the derivatives. Following 
Vinokur (1974), the strong conservation law form can be recovered by first dividing by the Jacobian then 
adding and subtracting like terms. For example, the term becomes 


E^, 

J 





Doing this for all the teims, and rearranging, results in 



(2.15) 


The last three terms, in braces, are called the metric invariant terms. By using the expressions for the metric 
coefficients, given by equations (2.13), one can show that the metric invariants are identically zero. In two 
dimensions, this is also true when derivatives are approximated by the finite difference formulas of Section 
5.O. 3 With the metric invariant terms eliminated, no metrics or flow variables appear as coefficients, and the 
strong conservation law form of the governing equations has been recovered. 


Equation (2.15) can be rewritten as 


dQ dE dF _ v 

d-t dl dr, d£ dr, 

where 



E = (££* + F£y + Q£ r ) 


(2.16) 


3 This is not necessarily true in three dimensions, however. 
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F = - j (Er) x + Ff) y + Qrji) 

Ev = j<£vt x + F v ty) 

F[/ = -j(E l />; JC + F I />/ > ,) 

Using equations (2.2a) through (2.2e) these can be expanded as 

Q = -jl> pw pv £ r ] r 

pw£ x + pv£j, + pf, 

(pu 2 + P)^ + pUV^y + put' 

puv£ x + (pv 2 + p)£ y 4- pv£ f 
(£ r + pK* + (£7- + + e t Zt 

purf x + pvrjy + prj t 
(pu 2 + p)t) x + puvtjy + purit 
puvr\ x + (pv 2 + p)r/ y + pvr, t 
(E T + p)ur\ x + (E T + p)vyj y + E T f] t 



where 



1 

1 

T xx%x 

F T xy%y 

V- 

J 

Re r 

*xy£x 

+ T yyly 




Pxtx 

+ Py£y 




- 

0 


1 

1 


F t jc yVy 

v = 

J 

Re r 

r xyV x 

+ fyyHy 




Px*lx 

+ PyVy 

Px 

= ur xx + 

v ?xy~ 

1 

Pr r qx 

Py 

= ur xy + 

VTyy~ 

1 

Pr r q y 


(2.17a) 


(2.17b) 


(2.17c) 


(2. 17d) 


(2.17e) 


In the viscous terms, the shear stresses and heat fluxes are defined exactly as in equations (2.3), except 
the derivatives in the Cartesian coordinate directions must be evaluated using the chain rule. For example, 

du _ du , du 
dx d£ x drj ^ x 

A A A A 

Note that F and F K have exactly the same form as E and Ev, but with f replaced by rj. 
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3.0 TIME DIFFERENCING 


The governing equations are solved by marching in time from some known set of initial conditions using 
a finite difference technique. The time differencing scheme currently used in Proteus is the generalized 
scheme of Beam and Warming (1978). The time derivative term in equation (2.16) is written as 


dQ AQ” gi d(AQ n ) i dQ n fl 2 AQ* ~ 1 
d-c At l + 62 dt 1 4- 6 2 dt 1 + 6 2 At 

+ c> [( 0 ]- i _02 ) AT,(AT)2 ] 


or, 


fr, BjA T a(AQ") At <?Q» 

y i + e 2 3 t + i + e 2 d-c 


AQ" ■ 1 + o (e, - \ - 0 2 )(At) 2 , (At) 3 


(3.1) 


where AQ" = Q" + 1 — Q". The superscripts n and n + 1 denote the known and unknown time levels, re- 
spectively. 

The parameters 6 1 and determine the type of time differencing scheme used. Some of the methods 
available with the above formula are given in the following table. 


0. 

*2 

Method 

Truncation Error 

0 

0 

Euler explicit 

0{AtY 

0 

- 1/2 

Leapfrog explicit 

0{ At) 3 

1 

0 

Euler implicit 

O(At) 2 

1/2 

0 

Trapezoidal implicit 

O(At) 3 

1 

1/2 

3-point backward implicit 

0(At) j 


Note that even though the generalized time differencing formula includes explicit methods, the Proteus code 
assumes an implicit method is being used. Note also that the truncation error listed in the table is the error 

A 

in the expression for AQ". The overall numerical method used in modelling the differential equations re- 

A 

quires AQ^/Ar, so the order of the overall method is this truncation error divided by At. 


A A 

Solving equation (2.16) for dQ/dr and substituting the result into equation (3.1) for c^AQ^/dT and 
dQ'/dr yields 

AO"= 6 ' AZ ( g (AF”) \ At / dt dt ) 

V 1+^y dt dr, J l + e 2 \ dt d v J 

0, At / d(AE y) d(AF v n ) \ At ( dE y n dE v n \ 

+ 1 +0 2 l dt + dr, ) + l+0 2 l dt + dr, I 

+ AQ” ~ 1 + o[(d } - ± - 0 2 )(At) 2 , (At) 3 J (3.2) 
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3UKI 



4.0 LINEARIZATION PROCEDURE 


4.1 INVISCID TERMS 


AAA A 

Equation (3.2) is nonlinear, since, for example, AE" = E n + l — E" and the unknown E n + 1 is a nonlinear 
function of the dependent variables and of the metric coefficients resulting from the generalized grid trans- 
formation. The equations must therefore be linearized to be solved by the finite difference procedure used 
in Proteus . This is done by expanding each nonlinear expression in a Taylor series in time about the known 
time level n. Letting G represent any nonlinear expression, 


G n + X = + At + 0(At) 2 (4.1) 

where 

dG _ dG dp dG djpu) gp d(pv) gG d£ r 

dr dp dr d(pu) dr d(pv) dr dEp dr 


Note that for linearization purposes only the metric scale coefficients have been assumed to be locally inde- 
pendent of time. Note also that for this linearization procedure to be second-order accurate, dG/dr (and 
therefore dpjdr , 3(pw)/<3r, etc.) need only be first-order accurate. Using forward differences, then, so that 



n + 1 

P 


n 

- P 


At 


+ O(At) 


A n 

At 


+ O(At) 


etc., equation (4.1) becomes 

G * * 1 - + ( f ) V + ( )" iM " + ( ag) )'**■ * ( M )"^ + (42) 


As an example the d(puv£ y )ld£ term from the x-momentum equation (part of the second element of 

A 

dEjd^) will be used. The nonlinear part of this term is ( puv) n + l . Rewriting this in terms of the dependent 
variables, 


(puv) 


n + 1 


jpuKpv) j 


Using equation (4.2), this is linearized as 

(puvf + 1 = {puvf - ( uv) n (p n + 1 - p n ) + v”[(pu) n + 1 - (pu)"] + «"[(pv) n + 1 - (pv)”] + O(At) 2 
which can be rewritten as 


A(puvf — — {uv) n Ap n -f v n A{pu) n 4- u n A{pv) n 4- 0(Ar) 1 


This linearization procedure, when applied to the entire AE" term in the vector equation (3.2), can be 
written as 
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n 


AE 


n 


^7T ) AQ" + 0 { At) 2 

5Q J 


(4.3) 


A A 

where (dE/dQ) n is a Jacobian coefficient matrix (not to be confused with the Jacobian J of the generalized 
grid transformation.) A similar equation can be written for AF". 


Each term in each element of E and F, given by equations (2.17b) and (2.17c), is linearized using the 

A A A A 

above procedure to generate the elements of the Jacobian coefficient matrices 5E/3Q and 5F/5Q. (Note 
that 3E /3Q = Jd E/5Q.) When this is done dEjdQ can be written as 


A 

dE 

A 

0Q 


ii 

dp . f 


it+A + u *x + 


dp 


d(pu) 


dp 

r Zjl 


0(pv) 


J?- z 

dE-r 


8p_ 

dp 


Zy ~ v f\ 




dp - 

v ^ x + ~d(JU) 

dp 


it+A + v Zy + 


JE-t 

d(pv) ^ 


dp 


(4.4) 


L 


fex +/i 


d{pu) 


Ky + f'S) 




where = u£ x + v£ y and f 2 = (£V + p)fp . The Jacobian matrix dFjdQ has the same form as dE/dQ, but 
with £ replaced by rj. 

The linearized pressure terms have deliberately been left in terms of dp! dp, dpjd(pu), etc. The ex- 
pressions to be used for these derivatives depend on the equation of state. Those currently built into the 
Proteu ; code, for a perfect gas, are presented in Section 4.3. 

4.2 VISCOUS TERMS 


The nonlinear viscous terms in equation (3.2), involving AEy and AFy, must also be linearized. To do 

A A 

this, the elements of Ey and Fy, given in equations (2.17d) and (2.17e), must first be rewritten in terms of 
the dependent variables, and with derivatives in the Cartesian directions transformed to derivatives in the 
computational directions using the chain rule. When the resulting expressions are substituted into equation 
(3.2), mixed second derivatives appear as well as second derivatives in a single coordinate direction. The 
mixed, or cross, derivative terms would lead to considerable complications in the implicit numerical solution 
algorithm if they were linearized using the procedure presented in Section 4.1. The two types of second 

A A 

derivatives are thus treated differently, and Ey and Fy are written as 


AAA 

Ek = E^, + 




(4.5) 


A A A A 

where Eyj and Fy, only contain derivatives in the £ and rj directions, respectively, and Ey 2 and Fy 2 contain 

A A 

derivatives in the other direction. The fully expanded expressions for Ek p Ek 2 , etc., are fairly long, and 
therefore are presented in Appendix A. 
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4.2.1 Non-Cross Derivatives 


A 

Examination of the elements of Eki in equations (A. 2a) through (A. 2c), and (A.2e), shows that every 
term has the form/^, where g is a function of the dependent variables, and/is a function of p, k, and/or 
the metric coefficients. Expanding in a Taylor series about time level n gives 

</Sj)" + ' - </*;)" + [ 1 At + 0(At) 2 

For linearization purposes only , we will assume /is locally independent of time. We can thus write 

<Ai>* - 1 - (/«()" +/" [ £ ] At + <K 4t) 2 

where 

dg = dg_dp dg d(pu) 
dr dp d-r d(pu) dr 


Therefore 

+ 1 = (/g t )n +/ ” ~k [' ^ Ap + a^4 A(pu) + " ] + 0(At)2 

A 

As with the in viscid terms, the linearization procedure for the entire AE£, viscous term in equation (3.2) can 
be written as 



AQ n + 0 ( At) 2 


(4.6) 


, A A A 

A similar equation may be written for AF^. The Jacobian coefficient matrix dEyJdQ is 



where 


(4.7) 


a xx = ( 2 M + t)i 2 x + piy 
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<*yy = PS X + ( 2 M + ))Zy 


a xy — 0* "F 




Like the pressure terms discussed earlier, the form of the temperature terms will depend on the equation 
of state being used. Those currently built into the Proteus code, for_a perfect gas, are presented in Section 
4.3. 

Note that in equation (4.6) the derivatives appearing in the Jacobian coefficient matrix 5E Kl /dQ are also 
to be applied to the AQ" appearing outside the parentheses. For example, the element in the second row 

and second column of dEvJdQ, which corresponds to the A (pu) term in the x-momentum equation, is 
a XI d{llp)ld£. For this term, the notation used in equation (4.6) means 

(- 5 -) 

\ / 22 

n J_( A(pt//7) n 
“ axx 3M P n 



The Jacobian coefficient matrix for the remaining non-cross derivative viscous terms, dFvJdQ, has the 

A A 

same form as <3E vJdQ, but with ? replaced by rj. 

4.2.2 Cross Derivatives 

As stated earlier, linearizing the cross derivative viscous terms in the same way as the remaining terms 
is very complicated wi thin the framework of the implicit numerical solution algorithm used in Proteus. 
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They are therefore simply lagged (i.e., evaluated at the known time level n and treated as source terms.) 
As noted by Beam and Warming (1978), this does not lead to a formal accuracy loss since 


AE y 2 = AE y~ 1 + O(At) 2 
AFt- = AF^r 1 + O(At) 2 


(4.8) 


4.3 EQUATION OF STATE 


The expressions to be used for dpjdp, dT\dp, etc., which arise from the linearization procedure, depend 
on the equation of state. The equation currently built into Proteus is for perfect gases, and can be written 
as 

P = { y-l )[E T -j-p(u 2 + v 2 )\ (4.9) 

or, in terms of temperature, as 

r . i [-^-|( u 2 + v ! )] ( 4 io ) 


With this equation of state, then, the appropriate derivatives are 


dp_ 

dp 



d{pu) 


= - (y - i )« 


<?P 

d(pv) 


= - (y - l)v 


ST 

dp 


dp_ 

dE T 


= y - 1 


1 

Cy 

~ J 

2 ~ 
L p 

to 

dT 

U 

d(pu) 

CyP 

dT 

V 

d{pv) 

C v p 

dT _ 

j_ 


dE r 


c,p 


(4.11a) 
(4.11b) 
(4.11c) 
(4. lid) 
(4.12a) 

(4.12b) 
(4.12c) 
(4.1 2d) 


If constant stagnation enthalpy is assumed, as discussed in Section 2.2, the appropriate equation of state 
is 

P = 1 ^-p[Ar-y(“ 2 + v 2 )] (4-13) 

and the temperature becomes 

T = -^[ h T-j(“ 2 + v 2 )] ( 4 - 14 ) 
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(4.15a) 


With these equations, the derivatives of p and T with respect to the dependent variables are 


dp v - 1 

d{pu) ~ y “ 


dp _ _ y - 1 

<?(pv) ~ y 


dT 

dp 


1 

CpP 


, 2 . 2 , 

(w +v ) 


dT u 

d{pu) c p P 


6T _ v 
5(pv) c pP 


4.4 LINEARIZED GOVERNING EQUATION 

The linearized form of equation (3.2) can now be written as 



(l+0 3 )Az f dE^ dFy 2 ^ 0 3 At f dEy^ ^ dF^, ^ 

+ 1+02 l d% + dt 1 J 1 + 0 2 y + dr 1 J 

+ AQ* - ' + 0^0, - - 0 2 )(At) 2 , (0 3 - «i)(At) 2 , (At) 3 J 


(4.15b) 

(4.15c) 

(4.16a) 

(4.16b) 

(4.16c) 


(4.17) 


There are a couple of things that should be mentioned about this equation. First, this equation is in 
so-called "delta" form. We will actually be solving this equation for AQ" and recovering Q” + 1 from 

Q" + 1 = AQ" + O' And second, in the coefficients of the cross derivative viscous terms the time differencing 
parameter 0i has been replaced by 0 3 . For second-order time differencing (i.e., if 0i = 0 2 + 1/2), 0 3 should 
be set equal to 0 3 . For first-order time differencing, however, 0 3 can be set equal to zero without losing 
accuracy. 
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5.0 SPACE DIFFERENCING 


To solve equation (4.17) an evenly spaced grid is defined in the computational (£, rj) coordinate system. 
Spatial derivatives are then approximated by finite difference formulas. First derivatives in the £ direction 
are approximated using the following second-order central difference formula. 



hkr 


fi+\J 


(5.1) 


The subscripts / and j represent grid point indices in the £ and r\ directions. The computational grid spacing 
is constant, and equal to 1/(M — 1), where N\ is the number of grid points in the f direction. A similar 
formula is used for first derivatives in the rj direction. 


The non-cross derivative viscous terms in the £ direction in equation (4.17) all have the form 

[/^T 


where Q represents one of the elements of Q. Using central differences this is approximated by 


_d_ 


\f-§£ C?A0)1 - S^fS^AQ^j 


££ ( fi + 1 /2,Ate A 0i + 1/2,7 — fi- 1 /2,A^0| - 1 /2,y'} 
-h- W+ 1/2 J(g*Q)i + 1,7 - feA0,,y] 


(A 


-fi- mjlte&Qkj- (sm~ 1 , 7 ]} 


2 m 


2(A f) J 


T «Ay + l,y)CfeA0 /+ u - feA0, y .] 
- <4y + ft- - teA0/_ !,;]} 

•{0/-w + 4;K?Afi)/-u 
~ (A- 1.; + 2 Ay + A+ i,;feA0i,7 

+ (Ay + fi+ l,y)O?A0 l + i,y} 


A similar formula is used for second derivatives in the r\ direction. 

Cross derivative viscous terms are evaluated using the following central difference formula. 


(5.2) 
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x 7 u 

- 2^ U+l ,/^r? ^)j + 1 J “ fi - 1 ,/^fj -S'); - 1 ,y] 

= 4A?A>7 1 >^ i+ W + 1 “&+W- l) 

fi — 1 ,yt?t — l,y + 1 — Si — 1 ,j — l)3 

Note that this formula is only needed for the source terms, since the viscous cross derivative terms are lagged 
one time level. 


When first derivatives are needed normal to a computational boundary, such as for Neumann boundary 
conditions, either first- or second-order one-sided differencing is used. The first-order formula at the = 0 
boundary is 


and at the £ = 1 boundary, 



1 


( flj-fxj ) 



</aW “ fv, - ij 


(5.4) 


(5.5) 


The second-order formula at the £ = 0 boundary is 

(df) ^ 3 f '>j + 4f2 ’j ~ /•/ 

x 7 ij 

and at the £ = 1 boundary, 

2A£~ ( /^i - W ~ 4 /v. - w + 3 /^i./ 

Similar formulas are used at the >7 = 0 and rj = 1 boundaries. 



(5.6) 


(5.7) 
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6.0 BOUNDARY CONDITIONS 


Choosing boundary conditions is perhaps the most important step in solving a flow problem with 
Proteus . Since the equations being solved at interior points are the same for every problem, the boundary 
conditions are what determines the final flow field for steady flows. 

With the difference formulas presented in Section 5.0, boundary conditions are required at each 
computational boundary, where N eq is the number of equations being solved. Note, however, that this is 
a numerical requirement, not a mathematical one. For example, for one-dimensional Euler flow N eq = 3. 
However, characteristic theory shows that, mathematically, only two conditions may be specified at a sub- 
sonic inflow boundary, and only one at a subsonic outflow boundary (Pulliain, 1986a). Some sort of ex- 
trapolation is typically used for the additional numerical boundary conditions. 

A variety of boundary conditions are built into the Proteus code, including: (1) specified values and/or 
gradients of Cartesian velocities u and v, normal and tangential velocities V n and V u pressure />, temperature 
T, and density p\ (2) specified values of total pressure pr, total temperature 7> , and flow angle; (3) linear 
extrapolation; and (4) spatial periodicity. Another useful boundary condition is a "no change from initial 
condition" option for u , v, p, 7, p, p T , and/or 7>. Provision is also made for user-written boundary condi- 
tions. The boundary conditions may be steady, unsteady, or time-periodic. The exact combination of 
boundary conditions to be used will depend on the problem being run. 

The boundary conditions in Proteus are treated implicitly. They may be viewed simply as additional 
equations to be solved by the ADI solution algorithm. And, in general, they involve nonlinear functions 
of the dependent variables. They must therefore be linearized using the procedure described in Section 4.0. 
The following sections describe this linearization for the general types of boundary conditions currently built 
into Proteus. 

6.1 NO CHANGE FROM INITIAL CONDITIONS, Aje=0 

This boundary condition simply sets the boundary value of the function g equal to its initial condition 
value. It can be written as 


A^ = g” + , -^ = 0 (6.1) 

A 

In general, g can be a nonlinear combination of the dependent variables Q. Linearizing g using the proce- 
dure described in Section 4.0, we get 


V dQ 


AQ” + O(A-r) 2 


( 6 . 2 ) 


Neglecting the 0(At) 2 linearization error, the linearized form of equation (6.1) can thus be written as 

n 

\ AQ” = 0 (6.3) 


dg 


\ dQ / 

6.2 SPECIFIED FUNCTION. e = f 

A specified function at a boundary can be written simply as 


f + '-f 


(6.4) 
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where g is the function being specified and /is the value being specified. Note that /can vary along the 
boundary, and can be time-dependent. Using equation (6.2) and neglecting the linearization error, the 
linearized boundary condition becomes 

^-^AQ n =/-/ (6.5) 


6.3 SPECIFIED COORDIN ATE DIRECTION GRADIENT. dzld<t> = f 


A specified gradient of a function in a coordinate direction can be written as 



»/ 


(6.6) 


where g is the function whose gradient is being specified, / is the specified value, and <p is the coordinate 
direction £ or y. Note that / can vary along the boundary, and can be time -dependent. 


The linearized form of g is given by equation (6.2). The linearized form of equation (6.6) can thus be 
written as 


Sg 

d<j> 



dg_ 

A 

dQ 



=/+ 0( At ) 2 


(67) 


Replacing differential operators with difference operators and neglecting the linearization error, the 
linearized boundary condition can be written as 


<5 


4> 




=f-6*g n 


(6.8) 


where <5* represents the one-sided difference operator to be used at the boundary. Options are available in 
Proteu r to use either first-order two-point or second-order three-point differencing. 


Note that this boundary condition is a specified value of the derivative with respect to the computational 
coordinate, not with respect to the physical distance in the direction of the computational coordinate. 
Following Kom and Kora (1968), and using the properties of the generalized coordinate transformation, it 
can be shown that for the £ direction the two derivatives are related by 


dg 


ds 


( 



dg_ 

dS 


Similarly, for the t] direction, 


dg _ J dg 

ds - “ /7T17 e ” 


If the value /= 0, of course, the two derivatives are equivalent. 

6.4 SPECIFIED NORMAL DIRECTION GRADIENT, Vg. n =f 

A specified gradient of a function normal to the boundary can be written as 

Vg” + 1 . n = / (6.9) 
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where g is the function whose gradient is being specified, /is the specified value, and n represents the unit 
vector normal to the boundary. Note that / can vary along the boundary', and can be time-dependent. 


For illustrative purposes, assume we are specifying a gradient normal to a constant £ boundary. Then 


n — ■ 


1L 

|V*| 


= ~^r z + 


1 . 1 r T 


m ** ‘ ‘ m 


where 


m 




Equation (6.9) can then be written as 


-^<sZ +, k+*; + %)“/ 


Using the chain rule to expand + 1 and + \ 


n + 1 n -hit . n + 1 

& £x + &, *lx 

n 4- 1 * + 1 t , n + 1 

Sy =•*{ ^ + ^ *1y 


Substituting into equation (6.10) and rearranging, 


is + 'ttl+ Zy) + £ + + W = rnf 


(6.10) 


Solving for g % + 1 , 


(If-) -4— + <61l) 

Now, in order to incorporate this equation into the ADI solution procedure used in Proteus , the dgjdrj term 
in equation (6.1 1) is lagged one level, and evaluated at time level rt instead of n 4- 1. Strictly speaking, this 
introduces an O(At) error into the solution. In practice, however, the actual error will depend on the degree 
of nonorthogonality of the coordinates near the boundary. For orthogonal coordinates no error is intro- 
duced. 


Using equation (6.2), and introducing difference operators and neglecting the linearization error, we can 
now write the linearized boundary condition as 


6 


f 


dg 

A 

dQ 




“T (ZxVx + Zy’1y)t v g n ~ S ^g n 
m 


(6.12a) 


where <5^ represents the one-sided difference operator to be used at the boundary. Options are available in 
Proteus to use either first-order two-point or second-order three -point differencing. 

Note that the unit vector n in equation (6.9) is in the direction of increasing £. Therefore, a positive 
value for /in equation (6.12a) indicates a flux in the direction of increasing £. Thus, a positive /at £ = 0 
implies a flux into the computational domain, and a positive / at f = 1 implies a flux out of the computa- 
tional domain. 


Specifying a gradient normal to a constant »/ boundary is done in an exactly analogous manner. The 
resulting equation is 
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6 n ( ) A Q" ~~m T + lyty 6 #* ~ V" 
\ dQ ) J "i 


(6.12b) 


where 


m = V v x 2 + n/ 


A positive value for /in equation (6.12b) indicates a flux in the direction of increasing y\. Thus, a positive 
/ at r\ = 0 implies a flux into the computational domain, and a positive /at tj = 1 implies a flux out of the 
computational domain. 

6.5 LINEAR EXTRAPOLATION 

Linear extrapolation from the two adjacent interior points is also available as a boundary condition. 
At the £ = 0 boundary, where / = 1 , this can be written as 


n + 1 n 4- 1 . n + 1 n 

Si - 2&+1 +S i + 2 =° 


(6.13) 


Note that this is equivalent to setting (Pg/dt*), + 1 = 0. Using equation (6.2), we can write the linearized 
boundary condition as 


(6.14) 


-2- AQf-2 -if- 4Q? + ,+( Hr I 4Q," +2 --«’ + 2i" +1 - 2 r +2 

«Q/, V<JQ/ (+ , V «}/ (+2 


Analogous extrapolation boundary conditions can easily be wntten for the remaining boundaries. 
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7.0 SOLUTION PROCEDURE 


7.1 ADI ALGORITHM 


The governing equations, presented in linearized matrix foma as equation (4.17), are solved by an alter- 
nating direction implicit (ADI) method. The form of the ADI splitting is the same as used by Briley and 
McDonald (1977), and by Beam and Warming (1978). Although the split equations can be developed in 
more than one way, in this discussion the approximate factorization approach is used. 

Letting LHS(4.17) represent the left hand side of equation (4.17), we can write 


LHS(4.17) 



AQ" (7.1) 


where I represents the identity matrix. Note that in this equation, using the djd £ term as an example, the 
notation used is meant to imply 


_d_ 

d( 




The term in braces in equation (7.1) can be factored to give 


LHS(4.17) 


1 + 


#iAt q 

i +e 2 di 


dE dEV, \ 

A A I 

dQ dQ ) 


#iAt a 

1+02 0*7 


dr 

A 

dQ 


d¥ Vi 

A 

dQ 




(7.2) 


The last term represents the splitting error. Note that, since AQ" = 0( At), this term can be neglected 
without affecting the overall time accuracy of the algorithm, even when second-order time differencing is 
used. 
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Equation (4.17) can thus be rewritten in spatially factored form, and, neglecting the temporal truncation 
and splitting error terms, becomes 


1 + 


A 

dE 


dE V] 


0,At <5 


\+e 2 ^ I 3 q d Q 


i + 


AjAx d 
1 + 3*1 


dT 5F V, 


At 


1 + 0 2 \ dri 


dE + dE_ } + At 


0,At ( 3E V 


1 @2 

n- 1 


dE V] 

dZ 


dE V] 

~dv 


+ 


AQ" = 

dQ dQ 

(1 + 0 3 )At ( dE Vi 0F 


A \" 
V , \ 


1 + 0 , 


dZ 


dr, 


1 + 0 , 




■ + • 


dE v . 


+ ■ 


1+0 


2 AQ" “ 1 


(7.3) 


Equation (7.3) can be split into the following two-sweep sequence. 
Sweep 1 direction) 


A, 

AQ + 


0] Ax d 
l+0 2 dZ 



0]At g 
l+0 2 dZ 



+ 


( A A \ n 

dE dF \ At 

di dr, J l+0 2 


dE Vi 

dZ 


(1 +0 3 )At 
1 + 0 2 


3Ey 2 d¥y 2 \ # 3 At 
d~Z + dr ] I ~ l + 0 2 




AQ 


n - 1 


(7.4a) 


Sweep 2 (n direction) 


AQ" + 


0j At 
1 + 0 2 


d_ 

dr, 


A 

_0F 

A 

dQ 



0]At d 
1 -h 02 3*1 



(7.4b) 


In the above equations, Q* represents an intermediate solution to the governing equations. 4 It should be 
noted that in Proteus, physical (i.e, n + 1 level) boundary conditions are used during the first ADI sweep. 
This introduces an <9(At) error in dQjds on the boundary for unsteady flows, but no error for steady flows. 
This point is discussed in detail by Briley and McDonald (1980). 


a a a 

4 The notation here is somewhat inconsistent The quantity AQ" = Q n + 1 — Q\ but AQ = Q Q\ not 

Qn+l _ Q* 
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Applying the spatial differencing formulas of Section 5.0 results in 


Sweep 1 (g direction) 


AQ f + 


0,At 


(1+0 2 )2A^ 

0, At 
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AQ, - 1 
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[u_ , + . aqI- , - (A- . + Vi +/ i+ ,)VAQ* + a + y; + ,)V+ iAQ* + ,] = 


(1 + 0 2 )2(A£) 

- t ^7 + s fT + ife^ (*&. + 6 ^ n 

(1 +0 3 )At 
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tl -t-t7 3/) iiT u A c A 0 3 At , £ A , £ A ^„_1 , e 2 A ^, 
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(7.5a) 


Sweep 2 (yj direction) 


a 0 T At 

AQ ' + (1 + 8 2 )2Arj 


dF 
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dQ 
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a q;+ , - 


;+i 


A 

dF 

A 

dQ 


AQ”-, 


y- 1 


0 1 At 


(1 + 0 2 )2(Aj,)" 


• [(£ - 1 + fj fgj- .AQ”- 1 - (j5_ 1 + 2£ + J5 + ,)' VaQ? + (fj + f j + ,)' V+ > A Qy + i] = 


AQ 


(7.5b) 


The subscripts / and j represent grid point indices in the £ and rj directions. For notational convenience, 
terms without an explicitly written i or j subscript are understood to be at i or y. In the viscous terms on 

A A 

the left hand side, / is the coefficient of d/d{ (or d/drj , depending on the sweep) in the dE v JdQ (or 

A A 

dFvJdQ) Jacobian coefficient matrix. Similarly, g is the term in the parentheses following djd £ (or d/Brj) 

A A A A 

in the dEyJdQ (or dFyJdQ) Jacobian coefficient matrix. Equations (7.5a) and (7.5b) represent the two- 
sweep alternating direction implicit (ADI) algorithm used to advance the solution from time level 
n to n + 1 . 


7.2 MATRIX INVERSION PROCEDURE 
7.2.1 Non-Periodic Boundary Conditions 


The complete set of algebraic equations for the first ADI sweep with non-periodic boundary conditions 
can be written in the following block matrix form. 5 


5 Although this discussion is written for the first ADI sweep, an exactly analogous procedure is followed for the 
second sweep. 
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These equations result from the application of equation (7.5a) for / = 2 to Ni — 1, with boundary conditions 

added at / = 1 and i = AT The parameter AQ* is the ^-element vector containing the unknown dependent 
variables; A, B, and C are the N eq x N eq coefficient submatrices at /— 1, U and / + 1, respectively; and S is 
the A^ 9 -element subvector containing the explicit source terms. Also, A', B\ and C' are the coefficient 
submatrices and S' the source term subvector for the boundary conditions. A variety of boundary condi- 
tions may be used. They are described briefly in Section 6.0, and in greater detail in Volumes 2 and 3. 


Note that the equations at the boundaries may contain coefficients at the boundary point and the two 
adjacent interior points. This occurs, for example, when extrapolation or second-order gradient boundary 
conditions are specified. As written, therefore, the coefficient matrix in equation (7.6) is not block 
tridiagonal. However, A{ can be eliminated by multiplying the second row of the matrix by A[ Cj 1 and 
subtracting from the first row. C'^ can be eliminated in a similar manner. Doing this, we define 

B, = B', - A'j CJ 'a 2 

C, -C'j -A'j Cj’Bj (7.7) 

Sj = S', - A', CJ % 

and 

a a’, = A V, ~ C'/v, A/v, 1 - 1 - 1 


®A r , — ® V, - C/Vj A Vl _ ,C^ _ , 
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The set of algebraic equations solved during the first ADI sweep can now be written as 



Since the coefficient matrix is now block tridiagonal, the equations can be solved using the block matrix 
version of the Thomas algorithm (e.g., see Anderson, Tannehill, and Pletcher, 1984). The procedure can 
be summarized as follows: 

1. Define Di = Bp 

2. Compute Ei = Dr L Ci and AQ1 = Dr ! Sp 

3. For / = 2 to N\, compute 
D, = B, — A t E,_, 

AQi =Dr 1 (S,.-A ( AQ;._ 1 ) 

(Actually, E, is only needed for i = 2 to N\ — 1). 

4. Then, set AQ^ = AQ ' Nl . 

5. Finally, for i = — 1 to 1, compute AQ, = AQ' — E,AQ 1 + 1 . 

A A 

In the Proteus code, in step 2 Ei and AQI are actually obtained by solving DjEj = Ci and DiAQi = Si 

A 

using LU decomposition of D. A similar procedure is used to compute E, and AQ' in step 3. 

7.2.2 Spatially Periodic Boundary Conditions 

In computational coordinates a spatially periodic boundary condition in the £ direction may be repres- 
ented as shown in Figure 7.1. 6 

6 As in Section 7.2.1 , this discussion is written for the first ADI sweep, but an exactly analogous procedure is followed 
for spatially periodic boundary conditions in the second sweep. 
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Figure 7.1 - Spatially periodic boundary condition. 

The grid points along the i — 1 and i = N\ lines are "similar" in the geometric sense, and have the same 
flow solution. Therefore, for a spatially periodic boundary condition in the £ direction, Qi = Q^,. 

To implement this boundary condition, an additional set of points is added at i=M+ 1, setting 

Q # , = Q 2 . This allows us to use central differencing in the <* direction at i — N\, computing the coeffi- 
cients in the same way as at the interior points. 
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The resulting set of algebraic equations will consist of M — 1 equations (for i— 2 to S\), with N x + 1 
unknowns. The block coefficient matrix thus has N\ — 1 rows and + 1 columns, as follows: 



A* 


r -1 

AQ t 
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A, 



A 2 B 2 C 2 
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A 3 ®3 C 3 

A * 

aq 3 


s 3 


A* 



A4 B4 C4 

aq 4 

m 


S 4 

■ 

■ 

m 

= 

■ 

■ 

• 

A * 


■ 

A V, - 2 ®iV, - 2 _ 2 
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aq;,-, 
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A AT t 

A * 

A Qtf, 


S .vi 

L J 

AQy 1 + , 




These equations result from the application of equation (7.5a) for / = 2 to N x . As in the previous section, 

parameter AQ* is the iV, 9 -element vector containing the unknown dependent variables; A, B, and C are the 
N eq x N eq coefficient submatrices at / — 1 , z, and z + 1, respectively; and S is the A^-element subvector con* 
taining the explicit source terms. 

A A A A 

Since Qi = and Q 2 = + 1 , equation (7.10) can be rewritten with N x — 1 unknowns as: 



A * 
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®2 C 2 A 2 

aq 2 
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AQv,-, 


S,v, - 1 

C AT, A /V, B /v, 

A * 

AQ lV , 


S ^i 

- 

_ 


- 


An efficient algorithm to solve this system can be derived that is similar to the Thomas algorithm for 
block tridiagonal systems. The procedure can be summarized as follows: 

1. Define D 2 - B 2 and F 2 = C*,. 

2. Compute E 2 = D 2 ] C 2 , G 2 = Dj *A 2 , and AQ 2 = D 2 1 S 2 . 
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3. For i = 3 to N\ - 1, compute 
D, = B; - A ; E, _ , 

E ( = D if'C, 

F/= — Fj-jE/-! 

G^-D-'AA-, 


AQ' = Df 1 (S ; — A,AQ' _ , ) 


4. Compute 

Gjv, _ i = ,(0^ _ ! — A^ _ ,0^ _ 2 ) 

EjV, - 1 = A tf, - E/v, - 2E.VJ _ 2 

' V 1 - 1 

D,v, = B,v, ~ X! 

( = 2 


AQ' S 



5. Then, set AQ^ = AQ'^ • 

A A A 

6. Compute AQ a ., _ , = AQV] - 1 - _ ]AQAr,. 

7. Finally, for i = N x — 2 to 2, compute AQ, = AQ' — E,AQ,+ — G AQ a ,- 


In the Proteus code, in step 2 E 2 , G 2 , and AQ 2 are actually obtained by solving D 2 E 2 = C 2 , D 2 G 2 = A 2 , 

a A 

and D 2 AQ 2 = S 2 using LU decomposition of D. A similar procedure is used to compute E„ G„ and AQ' 
in step 3, and G^-i and AQ'^i in step 4. 

7.3 UPDATING BOUNDARY VALLES 

7.3.1 Non-Periodic Boundary Conditions 

With the ADI algorithm described in Section 7.1, if gradient or extrapolation boundary conditions are 
used for the first sweep, the boundary values from the first sweep must be updated after the second sweep. 
This point is easiest to illustrate by looking at the following figure. 
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Figure 7.2 - Updating boundary values for non-periodic 
boundary conditions. 


In Figure 7.2, a 5 x 5 grid is shown in computational space. The triangles represent grid points at which 
the intermediate values Q* are computed during the first ADI sweep. These include the boundary points 

A 

at q = 0 and £ = 1. The circles represent grid points at which the final values Q^* are computed during 
the second ADI sweep, including the boundary points at rj = 0 and tj = 1. If gradient or extrapolation 
boundary conditions are used during the first sweep, so that the boundary values depend on the interior 
values, then the intermediate values at £ = 0 and £ = 1 must be updated after the second sweep to be con- 
sistent with the final values at the interior points. 

To do this, after the second sweep the boundary condition equations are rewritten and solved at the q 
boundaries. At the £ = 0 boundary, 

B',”AQ” + CfAQ” + AfAQ” = Sf (7.12) 

The subscripts refer to the value of i, the index in the £ direction. This equation is applied for j = 2 to 
A 2 — 1 in the rj direction. For notational convenience, however, the subscript j has been omitted. 



(7.13) 

(7.14) 


Proteus 2-D Analysis Description 


7.0 Solution Procedure 39 


(7.15) 


AQft = (B'vj )" ’ (sft - Cft AQft _ 2 - Aft AQft _ ,) 

Finally, note from Figure 7.2 that new comer point values are never computed in the solution algorithm. 
To make the comer values consistent with the rest of the flow field, in Proteus the comer values of density 
p and total energy £ r are arbitrarily defined by linearly extrapolating from the two adjacent points in both 
the £ and rj directions, and averaging the two results. The comer values of the velocities are updated by 
doing the same type of extrapolation. Instead of averaging, however, the extrapolated velocity whose ab- 
solute value is lower is used. This was done to maintain no-slip conditions at duct inlets and exits. 

7.3.2 Spatially Periodic Boundary Conditions 

Updating boundary values from the first sweep is complicated somewhat when spatially periodic 
boundary conditions are used. 
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Figure 7.3 - Updating boundary values for periodic 
boundary conditions in the £ direction only. 


The situation for a periodic boundary condition in the £ direction but not in the rj direction is shown 
in Figure 7.3. The triangles again represent grid points at which intermediate values are computed, and the 
circles represent grid points at which final values are computed. As can be seen from the figure, the inter- 
mediate values at f = 0 must be updated after the second sweep to be consistent with the final values at the 

A A 

interior points. This is easily done by setting Qi = Qvj for j = 1 to A 2 . 
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Figure 7.4 - Updating boundary values for periodic 
boundary conditions in the r\ direction only. 


The situation for a periodic boundary condition in the tj direction but not in the £ direction is shown 
in Figure 7.4. In this case, the intermediate values at f = 0 and at £ = 1 must be updated after the second 
sweep. To do this, the same procedure described in Section 7.3.1 for non-periodic boundary conditions is 

A A 

used, but for j = 2 to N 2 instead of N 2 - 1 • Then, for the lower comer values, Qi.i = Qi„v 2 and 

A A 

Q*l, 1 = Q^l.^2- 
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Figure 7.5 - Updating boundary values for periodic 
boundary conditions in both the £ and ij directions. 

And finally, the situation for periodic boundary conditions in both the £ and r\ directions is shown in 
Figure 7.5. Like the case with periodic boundary conditions only in the f direction, the intermediate values 

at £ = n must be updated after the second sweep. This is again done by setting Qj = Qa ( for y = 1 to N 2 . 
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8.0 ARTIFICIAL VISCOSITY 


With the numerical algorithm of Section 7.0, high frequency nonlinear instabilities can appear as the 
solution develops. For example, in high Reynolds number flows oscillations can result from the odd-even 
decoupling inherent in the use of second-order central differencing for the in viscid terms. In addition, 
physical phenomena such as shock waves can cause instabilities when they are captured by the finite dif- 
ference algorithm. Artificial viscosity, or smoothing, is normally added to the solution algorithm to suppress 
these high frequency instabilities. Two artificial viscosity models are currently available in the Proteus 
computer code - a constant coefficient model used by Steger (1978), and the nonlinear coefficient model 
of Jameson, Schmidt, and Turkel (1981). The implementation of these models in generalized 
nonorthogonal coordinates is described by Pulliam (1986b). 

8.1 CONSTANT COEFFICIENT ARTIFICIAL VISCOSITY 


The constant coefficient model uses a combination of explicit and implicit artificial viscosity. The 
standard explicit smoothing uses fourth-order differences, and damps the high frequency nonlinear insta- 
bilities. Second-order explicit smoothing, while not used by Steger or Pulliam, is also available in Proteus. 
It provides more smoothing than the fourth-order smoothing but introduces a larger error, and is therefore 
not used as often. The implicit smoothing is second order and is intended to extend the linear stability 
bound of the fourth-order explicit smoothing. 

The explicit artificial viscosity is implemented in the numerical algorithm by adding the following terms 
to the right hand side of equation (7.5a) (i.e., the source term for the first ADI sweep.) 

(2) a £ < 4) At 

(V{ AjQ + V^Q) - [(V { A^) 2 Q + (V„A„) 2 Q] (8.1) 

where and zf> are the second- and fourth-order explicit artificial viscosity coefficients. The symbols V 
and A are backward and forward first difference operators. Thus, 

= Qi+ i — Q; 

V ? AjQ* = Q ( + 1 -2Q l - + Q f _, 

(V,A ? ) 2 Q, = Q i + 2 - 4Q ( + , + 6Q, - 4Q, _ , + Q, _ 2 

Equivalent formulas are used for differences in the rj direction. 

A few details should be noted at this point. First, the sign in front of the artificial viscosity term being 
added to equation (7.5a) depends on the sign of the "i" term in the difference formula. For damping, that 
term must be negative when added to the right hand side of the equations (i.e., explicit artificial viscosity), 
and positive when added to the left hand side (i.e., implicit artificial viscosity.) See Anderson, Tannehill, 
and Pletcher (1984) for details. Second, the terms being added are differences only, and not finite difference 
approximations to derivatives. They are therefore not divided by A£, etc. Third, the variables being dif- 

A 

ferenced are Q, not Q. As noted by Pulliam (1986b), scaling the artificial viscosity terms by \ jj makes them 
consistent with the form of the remaining terms in the equations. Fourth, the terms are also scaled by At. 
This makes the steady state solution independent of the time step size (Pulliam, 1986b). And finally, note 
that the fourth-order difference formula cannot be used at grid points adjacent to boundaries. At these 
points, therefore, the appropriate fourth-order term in expression (8.1) is replaced by a second-order term. 
Thus, for points adjacent to the f = 0 and f = 1 boundaries, - e^ArQV^A^Q]// is replaced by 
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(8.2) 


+ 


•PAt 

J 


V^Q 


A similar expression is used at points adjacent to the rj — 0 and y = 1 boundaries. 

The implicit artificial viscosity is implemented by adding the following terms to the left hand side of the 
equations specified. 


- [V { Aj(/AQ*)] 


to equation (7.5a) 


£/A T 

“ 7 “ 


0 W yA Q”)] 


to equation (7.5b) 


(8.3) 


Note that the addition of the artificial viscosity terms, in effect, changes the original governing partial 
differential equations. At steady state, the difference equations with the artificial viscosity terms added ac- 
tually correspond to the following differential equations. 7 


dE dF 

dt dr, 


A 

dEy 

( 2 ) 

dFv ty 

dt 

dr, J 
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t E 

4 dVQ) 
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(4{) 



3*VQ) 

d{ 2 


+ (A*/) 2 


+ (An ) 4 


a 4 (VQ) 


aVQ) 

dr, 2 


The implicit terms do not appear, since they difference AQ, and in the steady form of the equations 

AQ = 0. The artificial viscosity terms do not represent anything physical. The coefficients should therefore 
be as small as possible, but still large enough to damp any instabilities. Although optimum values will vary 
from problem to problem, recommended levels are tf> = 0(1) and ti = 2s$p (Pulliam, 1986b). The recom- 
mended level for when used, is e?’ = 0(1). 


8.2 NONLINEAR COEFFICIENT ARTIFICIAL VISCOSITY 


The nonlinear coefficient artificial viscosity model is strictly explicit. Using the model as described by 
P ulliam (1986b), but in the current notation, the following terms are added to the right hand side of 
equation (7.5a). 




(- 7 ) + ( ~j ) J( e f > 4 {V - 4'VsVftJ 

(t) +(t). 

X 7 / + 1 V 7 J 


(4\q-4Va Q); 


(8.4) 


The difference operation ArV^A^Q is given by 

Af^jAfQ/ = Q; + 2 “ + 1 + — 1 

In the expression (8.4), \j/ is defined as 

^ 4 * \j/y ( 8 - 5 ) 


7 These equations represent the use of the constant coefficient artificial viscosity model presented in this section. The 
nonlinear coefficient model to be presented in Section 8.2 is more complicated, but the same principle applies. 
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where and \j/ y are spectral radii defined by® 


4> x = - 


U | + aj $ x 2 + Yy 


x Pv = - 


V\ + aj 


2 , 2 
y\x + Vy 


a>7 


Here U and V are the contravariant velocities without metric normalization, defined by 

U=£t+ £ x u + 
v== Vt + V 


( 8 . 6 ) 


(8.7) 


and a = yJyRT , the speed of sound. 

The parameters and c {4) are the second- and fourth-order artificial viscosity coefficients. Instead of 
being specified directly by the user, as they are in the constant coefficient model, in the nonlinear coefficient 
model they are a function of the pressure field. For the coefficients of the £ direction differences, 

(sj 2> )/= k 2 At maxfc* + i» ^i-i) (8.8a) 


(ef); = max[0, k 4 At - (e^ 2) )J 


(8.8b) 


where 


Aj- 1 ~ 2 Pi+Pi - l 
A+i + 2 Pi + Pi-\ 


(8.9) 


Similar formulas are used for the coefficients of the y\ direction differences. 

The parameter a is a pressure gradient scaling parameter that increases the amount of second-order 
smoothing relative to fourth-order smoothing near shock waves. The logic used tc compute e (4) switches 
off the fourth-order smoothing when the second-order smoothing term is large. 

The parameters k 2 and *c 4 are user- specified constants. Like the coefficients in the constant coefficient 
model, the optimum values will be problem-dependent, and are best chosen through experience. Cases have 
been run with values of k 2 ranging from from 0.01 for flows without shocks to 0.1 for flows with shocks, 
and k 4 ranging from 0.0002 for flows computed with spatially constant second-order time differencing to 
0.005 for flows computed with spatially varying first-order time differencing. Pulliam (1986b) gives 
k 2 — 0.25 and k 4 = 0.01 as typical values for an Euler analysis. 

Like the constant coefficient artificial viscosity model, the nonlinear coefficient model requires special 
formulas near boundaries. To apply (8.4) at / — 2, is needed at /= 1. It is defined as 

(4 2 >) 1 = K2ATmax(a 2 , <7j) 

With the above definition, applying (8.4) at / = 2 and / = Ai — 1 requires a at / = 1 and / = AT They are 
defined as 


8 It should be noted that the grid increments and A rj in these definitions do not appear in the corresponding 
formulas presented by Pulliam (1986b). This is because the grids used by Pulliam are constructed such that 
A£ = Arj = 1, while in Proteus A£ = 1 /(.Vj — 1) and A^ = I/OV 2 — 1). The definitions used here for \J/ X and \j/ y re- 
sult in an artificial viscosity level equivalent to that described by Pulliam. 
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- Pi + - 5 P2 + 2 P\ 

1 Pi + 4/>3 + 5/^ + 2p x 

- />,Vj _ 3 + 4/> V) _ 2 - 5/fy _ 1 + 2/fy 
<T,V ’ 1 “ P.V, _ 3 + 4PiV, - 2 + 5 ^.v, - 1 + 2 ^', 

And, finally, applying ( 8 . 4 ) at i = 2 and / = A', - 1 requires A { V { A { Q at i = 1 and / = JV, - 1- There are 
numerous formulas that could be used. The ones currently in the Proteus code are 

AjVjAjQ] = — Q 5 + 5Q 4 — 9Q 3 + 7Q 2 — 2 Q, 

A^V^AiQ^i _ ] = Q Ni _4 — 5Qjv, _ 3 + 9Q /Vj _ 2 — , + 2 Q^ 
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9.0 TURBULENCE MODELS 


As noted briefly in Section 2.0, for turbulent flow the Reynolds stress and turbulent heat flux terms are 
modeled using the Boussinesq approach. An effective viscosity is thus defined as fi = + fi tJ where \ii is 

the laminar, or molecular, viscosity coefficient, and fi t is the turbulent viscosity coefficient. Similarly, an 
effective second coefficient of viscosity is defined as X = 2/ 4- and an effective thermal conductivity coef- 
ficient is defined as k = k t -1- k t . 

The turbulent coefficients must be computed using a turbulence model appropriate for the flow being 
computed. In Proteus , turbulence is modeled using either a generalized version of the Baldwin and Lomax 
(1978) algebraic eddy viscosity model, or the Chien (1982) low Reynolds number k-s model. 

9.1 BALDWIN LOMAX MODEL 


For wall-bounded flows, (i.e., boundary layers), the Baldwin-Lomax turbulence model is a two-layer model, 
with 


(dinner for y n <y b 

(9.1) 

{(Pouter for y n >y b 

where y n is the normal distance from the wall, and y b is the smallest value of y n at which the values of fi t from 
the inner and outer region formulas are equal. For free turbulent flows (i.e., mixing layers, jets, and wakes), 
Mr = (Mr)<>u*r. In the inner region, in addition to the Baldwin-Lomax model, an alternate expression first 
presented by Spalding (1961), and later by Kleinstein (1967), is also available. 

In a simple boundary layer analysis, with only one solid surface, the procedure for computing /i t is rel- 
atively straightforward. In a general Navier- Stokes analysis, however, any or all of the boundaries may be 
solid surfaces. If both boundaries in a given coordinate direction are solid surfaces, the turbulence model 
is applied separately for each surface. An averaging procedure is used to combine the resulting two 
profiles into one. If neither boundary in a given direction is a solid surface, the formulation for free turbulent 
flows is used. In addition, values of fi t are computed separately for both the f and vj directions. This results 
in two complete turbulent viscosity fields. Another averaging procedure is then used to compute a single 
value of at each point in the flow. 9 

9.1.1 Outer Region 


The outer region turbulent viscosity at a given £ or r\ station is computed from 

(Pt) outer = KCcpPpKlebF wake^ e r (9*2) 

where K is the Clauser constant, taken as 0.0168, C cp is a constant taken as 1.6, and p is the static density. 
The parameter F wake is computed from 


9 This discussion is for the most general situation. When the flow is expected to be predominantly in one direction, 
input parameters in the Proteus code should be used to specify that direction. 
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F = 
1 wake 


I ymax F m 


c, 


wkVdiff f 


y'max 


for wall-bounded flows 
for free turbulent flows 


where C wk is a constant taken as 0.25, and 


v m= 


v 


(9.3) 


where V is the total velocity vector. 

The parameter F max in equation (9.3) is the maximum value of 


F(y n ) = 




for wall-bounded flows 
for free turbulent flows 


(9.4) 


and y max is the value ofy„ corresponding to F max . It has been found that for wall-bounded flows the function 
F(y n ) can have two peaks. As noted by Kirtley (1987), using the second peak as the location of F max yields 
the best results. 


For wall-bounded flows, y„ is the normal distance from the wall. For free turbulent flows, two values 
of F„ ax and y max are computed - one using the location of | V\ as an origin for y n , and one using the lo- 
cation of \ V\ . The origin giving the smaller value of y mcx is the one finally used for computing y n , F, 
and y max . 


In equation (9.4), |fi| is the magnitude of the total vorticity, defined for two-dimensional planar flow 
as 



dv _ du_ 
dx dy 


(9.5) 


The parameter A+ is the Van Driest damping constant, taken as 26.0. The coordinate y+ is defined as 




p^y n 

Fw 


^ r w p w Re r 

Re r = K, 


(9.6) 


where w, = Jr„lp»Re r is the friction velocity, t is the shear stress, and the subscript w indicates a wall value. 
In Proteus, -r„ is set equal to | | . 

The function F Kltt in equation (9.2) is the Klebanoff intermittency factor. For free turbulent flows, 
Fxitb = 1. For wall-bounded flows, 


FfCleb 


— (C Kieb) min 4" El Kleb)min^ 



^Klebht \ 

ymax J 


(9.7) 


This factor accounts for the experimentally observed fact that, as the free stream is approached, the fraction 
of time the flow is turbulent decreases. In equation (9.7), B and Cxm are constants taken as 5.5 and 0.3, 
respectively. {C Kleb ) min is a constant normally equal to 0.0. However, when using the Baldwin-Lomax model 
to generate initial turbulent viscosity values for the Chien k-t model (discussed in Section 9.2), (CxM)mn is 
set equal to 0. 1 . This yields a small positive value for p, in the free stream, and has been found to minimize 
starting problems with the k-t model. 
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9.1.2 Inner Region 


The inner region turbulent viscosity in the Baldwin-Lomax model is 

(H') inne r=pl 2 \n\Re r (9.8) 

where / is the mixing length, normally given by 

/= k j/j( i - e ~ y IA ) ( 9 - 9 ) 

and k is the Von Karman constant, taken as 0.4. 

A modified form of equation (9.9), proposed by Launder and Priddin (1973), may also be used. This 
formula is most useful for flows with steep negative gradients of shear stress normal to the wall, such as 
accelerated flows or flows with suction. Their modified formula for / is 

l=Ky„(l-e- y+{ ^ flA+ ) (9.10) 


where 


T + = 


T 




& 

W 


and n is a constant taken as 1.7. 


The inner region turbulent viscosity may also be computed using an alternate expression first presented 
by Spalding (1961), and later by Kleinstein (1967). In this model, 

(dinner = M/** - * 5 ^*" ~ 1 ~ ~ y ( xM+ ) 2 ] ( 9 - 1 0 

where 



T Pw^^r 


Again, in Proteus , t w is set equal to | £1 1 . 

9.1.3 Averaging Procedures for Multiple Boundaries 

As noted earlier, if both boundaries in a given coordinate direction are solid surfaces, the turbulence 
model equations are applied separately at each surface. It is assumed that the two inner regions do not 
overlap. The outer regions, of course, do overlap, and an averaging procedure is used to combine the two 
outer region profiles into one. For example, if the >7 = 0 and y\ = 1 boundaries are both solid surfaces, 10 
the two values of F wakt at a particular £ station aire combined using the following averaging formula: 

P (F W ake) 1 f\ (^wake) 2 fl /n 1 

Here {F wake )\ and (F wake ) 2 are the separate values computed for the rj — 0 and rj = 1 surfaces using equation 
(9.3). The parameters f x and f 2 are defined by 



10 An analogous procedure is used for solid surfaces in the £ direction. 
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/,= 


ggi 

Wi 


/r 


2D 7 

(y*h 


where n is a constant taken as 2.0, (y„)i and (Jr.): are the normal distances to the yj 0 and yj ^surfaces, 

respectively, and D t and D 2 are the normal' distances from the two >; surfaces to the location of m^. In 
addition, the y n }y m „ value needed in equation (9.7) for F K!eb is computed for both rj surfaces, and the mini- 
mum is used. These values of F wake and F K!cb are then used in equation (9.2) to compute 

The averaging procedure described above computes a single n, profile from the two profiles that are 
computed when both boundaries in a given coordinate direction are solid surfaces. We still must average 
the two values that result from computing (i, separately for both coordinate directions. 11 Following 
Goldberg and Chakravarthy (1987), this is done using the following formula: 


Mr = 


(Mr ly n ) 1 + (Mr tin) 2 tinhinh + tinhiFth 


[{Itinfl + ( 'tint]' 1 12 [tin)] + tin)l] 


2 - 11/2 


(9.13) 


Here (jr,), and (/r,) 2 are the separate values computed due to the presence of boundaries at £ = 0 and 
£ = 1, and at r\ = 0 and rj = 1, respectively. If there is only one solid surface in the £ direction, (y„)i is taken 
as the normal distance to that surface. If both £ boundaries are solid surfaces, (y„)i is taken as the normal 
distance to the closest one. If there are no solid surfaces in the £ direction, (y„)i is the normal distance to 

the location of either I V\ or I V\ , as described in Section 9.1.1. Analogous rules are used for (] y n ) 2 . 


9.1.4 Transition Model 

After fi, has been computed using the procedure described in the previous sections, a transition inter- 
mittency factor may be applied to simulate laminar-turbulent transition. The transition model is based on 
one given by Cebeci and Bradshaw (1984) for boundary layer analyses, and assumes that a geometric leading 
edge exists at either £ = 0 or y = 0. They report that the model is valid for adiabatic flows at Mach numbers 
less than 5. In this transition model. 


Mr = 


0 

yrrMr 


for x < x [r 
for x > x lr 


(9.14) 


where x is the distance from the leading edge, the subscript tr indicates a value at the start of the transition 
region, and y„ is a transition intermittency factor given by 


Vtr = 1 - ex P 


G(x - x„ 



(9.15) 


In equation (9.15), w, is the velocity at the edge of the boundary layer. The factor G is given by 

G = 8.33 x 10 -4 Re~ 134 

2 tr 

V 

where Re x „ = (u< xjv) rr Re n and v is the laminar kinematic viscosity at the edge of the boundary layer. 


11 As noted earlier, this discussion is for the most general situation. When the flow is expected to be predominantly 
in one direction, input parameters in the Proteus code should be used to specify that direction. 
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If we assume that, through the transition region, 14 ~ (24),, and v cs v fr , then equation (9.15) may be re- 
written as 


y t r= 1 “ ex P 


/ \ 2 

- 8.33 X 10 _4 ^ 66 ( - l) 


(9.16) 


To implement equation (9.16) in Proteus , we replace xjx tr with Re x /Re xtrl where Re x is defined as 


D 


Re x = 


Re r 


For flows predominantly in the direction, \V\ is the maximum total velocity magnitude at the current 
£ station, D is the distance from the point where j V\ = | V\ to the leading edge at £ = 0, and v is eval- 
uated at the point where I V I = I V I . An analogous definition of Re x is used for flows predominantly in 

the yj direction. 

9.1.5 Turbulent Values of / and k 


The turbulent second coefficient of viscosity is simply defined as 



(9.17) 


The turbulent thermal conductivity coefficient is defined using Reynolds analogy as 


*r = 


c p^i 

Pr t 


Pr r 


(9.18) 


where c p is the specific heat at constant pressure, and Pr t is the turbulent Prandtl number. In Proteus , the 
turbulent Prandtl number may be treated as constant, or as a variable using the following formula (Wassel 
and Catton, 1973): 



(9.19) 


Here C PrU C Pr2 , 0,3, and C P >4 are constants taken as 0.21, 5.25, 0.20, and 5.0, respectively, and Pr { = c p pLtjki 
is the laminar Prandtl number. 


9.2 CHIEN k-e TURBULENCE MODEL 
9.2.1 k-e Equations 

The low Reynolds number k-e formulation of K. Y. Chien (1982) was chosen because of its reasonable 
approximation of the near wall region and because of its numerical stability. Here k and e are the turbulent 
kinetic energy and the turbulent dissipation rate, respectively. 12 In addition, the Chien k-e turbulence model 
was frequently used in past Navier- Stokes computations with good results (Nichols, 1990, 1991; Patel, Rodi, 
and Scheuerer, 1985; Sahu, 1984.) The set of k-e equations are lagged in time and solved separately from 
the Navier-Stokes equations to allow for code modularity in turbulence modeling. In Cartesian coordinates, 
the two-dimensional planar equations for the Chien k-e model can be written using vector notation as 


12 It should be noted that in the Chien model, e is actually the isotropic portion of the turbulent dissipation rate. 
Throughout this manual, however, it is referred to as simply the turbulent dissipation rate. 
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(9.20) 


where 


and 


(9.21a) 

(9.21b) 

(9.21c) 
(9.2 Id) 

(9.21e) 


, Mr 

(9.22a) 

, Mr 

m £ = m + ^ 

(9.22b) 

C, = 1.35 

(9.22c) 

c,-c 2 ,( l-}e-* J ‘) 

(9.22d) 

a A = 1.0 

(9.22e) 

ct £ = 1.3 

(9.22f) 

II 

QO 

(9.22g) 

* - pkl 
K t ~ fit 

(9.22h) 

P “ = Re/'- 

(9.23a) 
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p '=((f ) 2+ ( 

' dv \ 2 (du 

K dy J j 3 \dx + 

\ 2 / v 2 

dv \ 1 du dv \ 

dy ) \ ^ dx ) 

(9.23b) 


P _ du i dv 


C9 73rt 


2 dx dy 



The turbulent viscosity is given by 

r * 2 

Pt = “ 


(9.24a) 


ii 

** 

1 

1 

n 

f 

) 

(9.24b) 


C M = 0.09 


(9.24c) 


C 3 = 0.0115 


(9.24d) 


Note that the vectors W, F, G, and S are used in most standard k-z formulations (with different con- 
stants), and the vector T is unique to the low Reynolds number formulation of Chien. The parameter y n 
is the minimum distance to the nearest solid surface, and y+ is computed from y n . The production of tur- 
bulent kinetic energy P k includes the full Boussinesq approximation for compressible flows. All of the 
above equations have been nondimensionalized using appropriate normalizing conditions. 
Nondimensionalization of mean flow properties is discussed in Section 2.1. The turbulent kinetic energy 
k and the turbulent dissipation rate z have been nondimensionalized by u? and p r i4hi r , respectively. 

Following the procedure of Section 2.3, the following generalized grid transformation is used to trans- 
form the k-z equations from physical (x,y, t) coordinates to computational (£, tj , t) coordinates. 

rj = rj(x,y) (9.25) 

T = t 


Applying the generalized grid transformation to equation (9.20) yields 

W T + F^ x + F n r, x + G^ y + G^y = S + T (9.26) 

Although the above equations can not be put into exact strong conservation law form, the procedure 
used to do so for the mean flow equations, described in Section 2.4, is nonetheless applied to equation 
(9.26). The result is 


where 


aw dF dG 

di dtf 


A A 


S + T 





A 



J 


<Z x puk 4- Zypvk 
t, x puz + Zypvz 


(9.27) 

(9.28a) 

(9.28b) 

(9.28c) 
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(9.28dJ 


A 



J 


1 

Re r 


ptfl + zfrt 


A 



J_ 

J 


1 

Re r 


^ki^x^x ^yVy)^rj 
f^ei^x^x 


A A A 

G = G r - G 


D ' 


'M 


Gr = 


J 


r] x puk + 

Y\ X pUZ + tfyPVE 


A 



A 1 

g a/=7 


s = 


1 

y 


1 1 

y Pe r 


Ski’ll + *h) k r, 

H e (*t 2 x + »?>K, 


1 

/?e r 


HkiZxlx + 
vMx*lx + Zyly)^ 


P k — /?e r />£ 


T = 


1 

y 


2 

/?e r 



2 



Re r 


2 


-V/i 


e 


(9.28e) 

(9.28f) 

(9.28g) 

(9.28h) 

(9.28i) 

'(9.28|) 

(9.28k) 


Note that in equation (9.28j), the term P* involves derivatives with respect to the Cartesian coordinate 
directions (see equations (9.23a-c.) These are evaluated using the chain rule. 

9.2.2 Linearization of the k-z Equations 


Solving equation (9.27) for dW/dr and substituting the result into the time differencing scheme of Beam 
and Warming (1978), given by equation (3.1), for d(AW ")|d^ and dW^/dr yields 


a 0, At 

AW = 1 


+ 


1 + 

At 


d(AF c ) | d(AF g ) [ ^(AF^) d(AG c ) < d(AG D ) | d(AG w ) < ^ ^ 


dt, 


A 

dF, 


di 


dZ 


dr, 


dr, 


dr, 


1 + 02 


- + ■ 


dF 


D 


dF 


M 


dZ dZ 


dZ 


A 

dG c 

drj 


- + 


3G 0 dG M 


dr, dr, 


+ S + T 


+ 7 Te; 4 ^" " ’ + °[( e ' - 2 " ,4T)J ] 

Equation (9.29) is then linearized using the procedure described in Section 4.0. Let 


A A 

d¥ c dF 

T = — B = 


D 


d W 


A > 

dW 


c= 


dG/ 


A ’ 

dW 


D -!Ra. m =JS- n= il 

U A » lVi A ’ 1 A 

dW dW dW 


(9.29) 


(9.30) 
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be the Jacobian coefficient matrices, where 


A = 


ZxU + ( y v 0 
0 c, x u + 


(9.31) 


B = 


JRe, + P ) 


JRe r P ) 


(9.32) 


C — 


Y\ x u + rjyV 

0 Y\ x u 


° 1 
+ vj 


(9.33) 


D = 


1 +iJ)(t) 


' n 


(9.34) 


A/ = 


-yr Jl£l 

£ M, 

P k £ 2 


* 2 ^ 

— C — — — /?p 

S* £ 2 ji f r 

-2Re r C 2 f 


(9.35) 


jV = 


2m 

py^Re r 

0 


0 

__Jp_ e ~y + l 2 
pyn Re r 


(9.36) 


The linearized form of equation (9.29) can now be written as 
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^v) _ w + a«*w ) _ Q(Zmw) . 

1 + 0 2 5f <3£ dt, Or, 


Mr d(AF„) a(AG jVf ) 
1 + $2 a >7 


AAA AAA 

aF r 5F n aF^ aG r aG^ aG ;V/ 


if- + if- *-sr tS+T + t^ 


4- AW” " 1 


> (fii-y-e 2 ) 


(At) 2 , (At) 3 


9.2.3 ADI Algorithm for k-s Equations 

Letting LHS(9.37) represent the left hand side of equation (9.37), we can write 


LHS(9.37) = jl + ^ A J- [-^■{A-B) + -~^(C-D)-(M + N)jl AW* 


(9.38) 


where I represents the identity matrix. Factoring the term in braces, and neglecting the temporal truncation 
and splitting errors, equation (9.37) becomes 


1 + -n^ [ i <" - + ^H 1 + [ i (C - fl> ] 


AW* = 


0,At d(AF w ) d(AG w ) 

T+T 2 di + dr, 


At ( d¥ c d¥ D d¥ M dG c dG D 8 G m a a\ 

+ T— H ^ 2 + xZ + » + * I t- 


di ac dt drj dr\ dt] 


(9.39) 


Using the procedure of Douglas and Gunn (1964), as written by Briley and McDonald (1977), equation 
(9.39) can be split into the following two-sweep sequence. 

Sweep 1 (a direction) 


0i At F ft I A * At 0(Ar »*) 

_J_[^- M . i) - ( , W + A,jt 4W -JJ- + 


.* At d(AF w ) d(AGtf) 

1 +S 2 si + *1 


ihL + ^L + i!f._i5£ + i^ + ^« + S + T) +t ^- 4W— (9.40a) 

d£ dQ d£ d?j dt] drj f 1 + t / 2 


Sweep 2 (tj direction) 


A A * 

AVV = AW 
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To approximate spatial partial derivatives, the spatial differencing formulas of Section 5.0 are used in 
equations (9.40a) and (9.40b). Following Nichols (1991), the spatial derivatives for the convective terms 
are approximated using first-order upwind differencing. A first-order backward difference approximation is 
used for the terms with positive eigenvalues, and a first-order forward difference approximation is used for 
the terms with negative eigenvalues. 


9.2.4 Matrix Inversion Procedure for k-e Equations 


Non-Periodic Boundary Conditions 


Explicit boundary conditions are used for ease of implementation and modification. The complete set 
of algebraic equations for the first ADI sweep with non-periodic boundary conditions can be written in the 
following block matrix form. 



In the above matrix, A„ B,, and C, are the 2x2 coefficient matrices in front of AW *_ u AW*, and 

A 

AW* +1 , respectively, and they should not be confused with the Jacobian coefficient matrices A, B , and C 
defined in equations (9.30). The above block tridiagonal coefficient matrix is solved using the Thomas al- 
gorithm discussed in Section 7.2.1. An analogous procedure is used for the second ADI sweep. 


Spatially Periodic Boundary Conditions 


A spatially periodic boundary' condition in the £ direction may be represented as shown in Figure 7.1. 

Following Section 7.2.2, an additional set of grid points is added at / = N\ + 1, setting W Nl + j = W 2 . This 
allows us to use central differencing in the £ direction at i = AV 


A A A A 

Since W, = Wat, and W 2 = YV v , + b equation (9.41) can be rewritten as: 
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To solve the above system, the algorithm described in Section 7.2.2 is used. An analogous procedure 
is used for the second ADI sweep. 

9.2.5 Updating Boundary Values for k-z Equations 

For easy modification and easy accommodation of complicated boundary' conditions for k and s, non- 
periodic boundary conditions are treated explicitly in the solver. After the k and t values at the interior 
points are advanced in time, the values at the boundaries are simply computed from the new interior values 
using the specified boundary conditions. 

Spatially periodic boundary conditions in either sweep direction are treated implicitly, as described in the 
previous section. For a periodic boundary condition in the f direction, the k and t values at / = 1 are easily 

updated by setting Wi = W^. An analogous procedure is used for periodic boundary conditions in the rj 
direction. 

9.2.6 Turbulent Values of A and k 

The turbulent second coefficient of viscosity ). t and the turbulent thermal conductivity coefficient k t are 
defined as described previously in Section 9.1.5. 
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APPENDIX A - EXPANSION OF VISCOUS TERMS 


In Section 4.2, the viscous terms in the governing equations are linearized. To do this, the elements of 

E v and FV, given in equations (2. 17d) and (2.17e) must first be rewritten in terms of the dependent variables, 
and with derivatives in the Cartesian directions transformed to derivatives in the computational directions 

A A 

using the chain rule. The non-cross derivative terms, involving E Vl and F^, are then linearized using Taylor 

A A 

series expansion. The cross derivative terms, involving E v 2 and Fk 2 , are simply lagged one time level. This 
Appendix presents the fully expanded viscous terms required in the linearization procedure. 


The viscous term Ek is given by equation (2.17d), which is repeated here. 



J 1_ 

J Re r 


0 

T xx < ’jc T r xy£y 
T xy%x T r yy^y 


Px%x + PyZy 


(A.l) 


where 

r xx = 2 P- U x + X ( u x + v y) 
lyy = 2 flV y + X(u x + V y ) 

T xy = v(Uy + v x ) 

Px = m xx + v ?xy--^r4x 
Py ^xy T VT yy p r fy 

4x=- kT x 

q y = - kT y 


The chain rule is used to transform derivatives in the Cartesian directions into derivatives in the com- 
putational directions, resulting in 

T xx = (2 M + + r\ x U v ) + + ri y v v ) 

Vyy (2/1 T 2)(CyV^ qyVfj) T ^ ( ^ T 

'xy = My u l + *ly u r, + + 1x v n ) 

Px = ( 2 M + -W^uuj + *1^) + 2(^uv ? + ti y uv v ) 

+ niQyVUt + rjyVUy + Z x Wr + r, x w v ) + — (^Tj + nj n ) 

py = (2m + 2)(^Wj + r, y w v ) + ;.(^v^ + rj x vu v ) 

+ niiyicu, + Vy uu n + £, x uv { + y,^) + — (ZyT^ + r, y T v ) 
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The above expressions for the t's and /?'s are next substituted into equation (A.l). The £ derivative 
terms become elements of E^, and the r\ derivative terms become elements of E y 2 . The resulting four ele- 


A 

ments of Ep-, (excluding the 1 Re, coefficient) are 

(E F[ ), = 0 (A.2a) 

(E Ki ) 2 = + >-Z x (Ax u > + + ^y^y u l + &*) ( A - 2b ) 

(E^) 3 = 2^)v l + + i y v { ) + rf x (ZyU> + £ x v ? ) (A. 2c) 

(E„,) 4 = 2 ii(guu s + ^Wj) + iSjLZxUUt + 

+ fJ.£ X {£yVUg + £ x w { ) + nSyityUUs + £ x z/v ? ) + yp- (£ 2 + £ 2 )T> (A. 2d) 

For linearization it is convenient to rewrite the last element as 

(E^) 4 = (2M 2 +X - [£ 2 (k\ + £j(v 2 ) s ] + (m + 2)£ x £ y (i/v) { 

+ y [&v 2 )« + (fx + #r { (A.2e) 

A A 

The elements of F n have exactly the same form as those of Ek„ but with £ replaced by rj. 

The four elements of E y 2 (again excluding the 1 \JRe, coefficient) are 

(E Ka ), = 0 (A. 3a) 

(E[/ 2 ) 2 = 2 tfxnxS + ttxiWr, + + vZyiVyUr, + *lx V r,) ( A - 3b ) 

(E K2 )3 = 2 + Kytojfy + VV + ^x(*ly u r, + »txV (A- 30 ) 

(E„ 2 ) 4 — 2fj.(^ x rj x uu r] + (.ytiyW^) + X^ x {y\ x uu^ + rtyuvj + ).^y(tj x vu^ + rjyW^) 

k 

+ fi£ x (riyVu n + r\ x w v ) + nZyfayUUq + ri x uv n ) + — (£ x >/ x + i y Vy)T n (A. 3d) 


A A 

The elements of F K2 have exactly the same form as those of E^ 2 , but with £ replaced by rj and r\ replaced 
by £. 
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APPENDIX B - AXISYMMETRIC AN ALYSIS 


The analysis used in Proteus for axisymmetric flow is essentially the same as for two-dimensional planar 
flow, described in the main body of this report. However, there are some additional terms in the 
axisymmetric equations that complicate things somewhat. For that reason, the axisymmetric analysis is 
described separately in this appendix. 

B.l GOVERNING EQUATIONS 

In cylindrical coordinates, the governing equations for axisymmetric flow, with swirl, can be written 
using vector notation as 


<XrQ) d(r E) <?(rF) d(r E v ) 

dt dx dr + dx 


d{r¥ v ) 

dr 


+ Hy 


(B.l) 


where 


Q = lp pu pv pw E r y 


(B.2a) 


E = 


pu 

pu 2 + p 

puv 


puw 

( [E r +p)u 


(B.2b) 


F = 


pv 

puv 

pv +p 


pvw 

(E T + p)v 


H = 


0 

0 

-p- pw 2 


pvw 

0 


(B.2c) 


(B.2d) 


Ei/ = 


Re r 


‘JOC 

T xr 

T x0 


m xx + VT xr + M x e - 7“ 


(B.2e) 
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0 


Yv = 


Re r 


r r6 


“?xr + VT rr + Wr r0-~p^ 




(B.2f) 


Hi/ = 


Re r 


0 

0 

_ T ee 

r r6 

0 


(B.2g) 


Equation (B.l) thus represents, in order, the continuity, x-momentum, r-momentum, 0-momentum (swirl), 
and energy equations, with dependent variables p, pu, pv, pw, and Er- Note that the additional terms in 
these axisymmetric equations destroy the strong conservation law form of the two-dimensional planar 
equations presented in Section 2.1. Unfortunately, the axisymmetric form of the equations cannot be put 
into strong conservation law form (Vinokur, 1974.) 

The shear stresses and heat fluxes are given by 


r ee = 2 m — + 



( d(rv) 

V Sr JJ 



k Sr )\ 

+{&+*( 

d(rv) \1 

^ )\ 

/ du . dv 
dr dx 

) 

dw 


II 

%i? 

1 

) 

. dT 



(B.3) 


<Jr = — k 


8T_ 

dr 


In these equations, x, r, and 0 represent the axial, radial, and circumferential directions, respectively; and 
u, v, and w represent the velocities in those directions. The remaining symbols are the same as those in the 
two-dimensional equations described in Section 2. 1 . 

For turbulent flow, p, and k represent effective coefficients. The turbulence model is described in 

Section 9.0. The only modification to the model for axisymmetric flow is the definition of | Cl | , the mag- 
nitude of the total vorticity. For axisymmetric flow, 
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Q. 


( dw _vv _ \ 2 / dw \ 2 / dv _ jht_ V 

V dr + r ) + \ d x ) + V dx dr ) 


When the generalized grid transformation of Section 2.3 (with y replaced by r), is applied to equation 
(B.l) the result is 

(r Q) t + (r Q)|£, + (r Q) TJ r tr + ( r E)^ x + (r E)^»j x + (r F) f + (r F) t) r lr + H 

- (r E V ) S Z X - (r E v ) vVx - (r F v ) g ( r - (r F v )^ r — Fly = 0 (B.4) 


Although this axisymmetric equation cannot be put into exact strong conservation law form, the pro- 
cedure used to do so for the two-dimensional equation, described in Section 2.4, is nonetheless applied to 
equation (B.4). The result is 


djrQ) d(r E) 
d-r + dl 


d( r F) 

dr, 


A 

+ H = 


djrEy) 

d{ 


d(rF y ) 

dr, 


A 

+ H y 


(B.5) 


where 


E = -Jr(E£ x + F£ r + <K,) 
F = y (Erj x + Frj r + Q»j f ) 


Ey — -J (Ey £ x + Fy £ r ) 
Ey = y (E y r, x + F y r, r ) 


Hy = 



Using equations (B.2a) through (B.2g) these can be expanded as 


Q = 


1 T 

— \_p pu pv pw E t ~] 


E 


1 

J 


pu£ x + pv£ r + pz t 

(pu 2 +P)Z X + puv£ r + pu£, 
puv£ x + (pv 2 +p)t r + pv£, 
puw£ x + pvw£ r + pw£ t 
(E T + p)ui x + (E r + p)v£ r + E T l l 


(B.6a) 


(B.6b) 
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(B.6c) 


A 

F 


pur\ x + pvr\ r + pri t 
{pu 2 + p)rj x + puvt lr + pltr l[ 
puvrix + ( pv 2 + p)y i r + pvr] t 
puwrj x + pvwrjr + P Wr h 
{E T + p)ur, x + (E T + p)vr! r + E T t] r 


H 


1 

J 


0 

0 

- p-pw 2 

pvw 

0 



I 

Re r 


0 

T xx^x + r xr%r 
T xr%x F r rr^r 
x x0^x + r re£r 

I Sxtx + PrSr 


0 



1 

Re r 


T xx t lx F T xr t lr 
z xr r lx F T rr r fr 

T x0Vx + T r0»/r 

fix*lx F fiplr 


H 


v — 


1 1 
J Re r 


0 

0 

— T ee 

T r0 

0 


where 

fix = UT xx F VT xr F WT X 0 p r tfx 
fir = m xr + VT rr + w ?r9 ~ 4r 

B.2 LINEARIZATION 

Solving equation (B.5) for cQjdx (assuming r is not a function of time) and substituting the 
time differencing scheme of Beam and Warming, given by equation (3.1), for d(AQ")/dr and 


(B.6d) 


(B.6e) 


(B.6f) 


(B.6g) 


(B.7) 


result into the 
dQ'/d r yields 
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aq"» ■( -3^+^4+ah ^+^ + h ; 


1 + 0 2 ' ^ dQ dr, 

0jAt ! / 3(rAEj/") d(rAFV") 


1+0, r 


<?£ 


dr, 


1+0, r 


d{ 


■ + ■ 


dr, 


+ ah/ + 


Ar 1 ( d(rE v n ) d(r F v n ) 


1 + 0, r 


dS 


+ 7 Te~ 2 4Q" ' 1 + o[(»i - -i- - e 2 )(A+, (A,) 3 ] 

This equation must be linearized using the procedure described in Section 4.0. 
B.2.1 Inviscid Terms 


A A 


For the inviscid terms the Jacobian coefficient matrix dE/dQ is 


dr, 


+ H. 


(B.8) 


dE 

A 

dQ 


dp 

dp 

dp 

dp 


l t 

lx “ u f\ 
lr ~ V /l 

~ W A 


lr +/j + ul x + 


dp 

d{pu) 


lr 

dp 

U lr lx 


dp 


dp 

v?1+ 5 7pZ) lr 

K, 

ap 


a(pv) 

? f +/l + v £r + 




a<pv) 


fr 




dp 

d(pw) 

lr+A 

dp 


lx 


lr 


dp 

dE r 

dp 

dE r 

0 


lr 


-/.(/.-•gr) r,-^ {,+/.(■ +-g^) 


(B.9) 


where /I = + v£, and = (£ r + p)lp. The Jacobian matrix dF/dQ has the same form as dE/dQ, but with 

^ replaced by r,. 

For the additional term H, the linearization procedure gives 


0 0 0 


0 0 


dH 

A 

dQ 


0 0 0 


0 0 


dp 2 dp 

5(p«) 

— vw 0 


dp d p 

d{pv) d(pw) W dE T 

w v 0 


0 0 0 


0 


0 


(B. 10) 


B.2.2 Viscous Terms 

A A 

To linearize the viscous terms, E Vv E v 2 , etc., must first be rewritten in terms of the dependent variables, 
and with derivatives in the cylindrical coordinate directions transformed to derivatives in the computational 
directions using the chain rule. The shear stress and heat flux terms, given by equations (B.3) and (B.7), 
become 
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T *x = (2 m + + Wr,) + T CW")s + »b( rv )»P 

T rr = 2 m(C,Vj + n r v,) + 2(^J + TJ X U V ) + -y \_$ r {rv)t + r,^)^ 
r ee = 2 flyr + + Vx u r,) + T [W")f + IrWr,] 

Txr = rtttfg + *>'“*) + ^ V ? + 


?re = + »?,%) -(“T 

/?* = (2m + A)(^ { + rtyaj + yr UM^k + *» A™)^ 

+ n(Z r vu> + r] r vu v + £ x Wj + njc%) + m(^hw ? + j;^) 

+ ^r(^ + ^) 

Pr = 2 M(£r w J + >7|- vv » r ) + 4“ KrK"){ + »?r v ( n ')^] 

+ M(^rWW{ + yj^ + £ X WV ? + »7 X %) + + tj.WW^) 

2 

+ 2(^V W| + t, x VU n ) - M-7- + -j^T (Zr T S + VrT'r,) 

The above expressions for the shear stress and heat flux terms are substituted into equations (B.6e) 
through (B.6g). As in the two-dimensional planar case, the cross derivative terms are separated from the 
non-cross derivative terms. In addition, for the axisymmetric case the non-derivative terms are included 
with the cross derivatives. 


The resulting five elements of E K , (excluding the \jJRe r coefficient) are 

(E Ki )! = 0 (B.lla) 

(E Ki ) 2 = 2/i^Wj + + yr ^(rv)^ + + £ x Vj) (B. 1 lb) 

(Ej/^3 = 2M^r v ^ + + yr J + ^x(^r u ^ + £jc v >) (®- * lc ) 

(Em,) 4 = ( Blld > 


(E„,) 5 = 2 ni&UUt + ^Wj) + + T W")*] + '^\}x vu l + T «,*("){] 

+ + i^vv { + ^ww,) + fi£ r (Z r uu ! ' + ^nv ? + ^vw ? ) + (£* + £?)T ? (B. 1 le) 

For linearization it is convenient to rewrite the last element as 

(Em,) 5 = (2M ^ - [&u\ + £(v\] + 0* + WxtMt + W 

+ y [^x( y2 + vy2 )| + ^ 2 (“ 2 + w 2 )?] + " f (£jc + (B. 1 If) 
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The elements of F Vl have exactly the same form as those of E^, but with £ replaced by rj. 


A 

The five elements of E y 2 (again excluding the 1 jJRe r coefficient) are 

(E^i=0 (B.12a) 

(Ey 2 ) 2 - 2^ x r, xUri + + T J + X r (*1r U r, + *)x v r,) (B. 12b) 

(E^ 2 )3 = 2 Xrtr v r, + Xr[*lx u r, + ~T »h( n ') J + X x {> b u r, + 1x v v ) (B. 12c) 

^ w 

(E^) 4 = X x *l x ™ v + n£rfr w r, - X r — (B. 12d) 


(Ej/ 2 )s = MZxlx*^ + frVS,) + + T + Hr\_Vx vu r, + T VW,] 

+ Xxil r™r, + Vx w v + + XXrVUr, + VxXr, + V*%) 

- ntr^J-+-fir(Zx r lx + Zi*lr) T r, (B.12e) 

The last element can be rewritten as 


( E ^) 5 = 2 M (^ x *i x uw n + trfrWj + k£ x (r ix uu ri + tj r uv v ) + Af r (»fjc v “», + ri r wj + Atj r —(Z x u + ^ r v) r n 
+ M-JjcOVH, + JJx% + r) x ww v ) + n£ r (r Ir UU rl + r,^ + i,^) 

2 


-Xr^ + -jt r (t x *lx+tr’1r)T n 


(B.12f) 


The elements of F Vl have exactly the same form as those of Ev 2 , but with £ replaced by >j and >; replaced 

by£. 


The five elements of Hk are 


(H„), = 0 


(H„) 2 = 0 

(Hf)3 = - 2m y - ;.(Z x Ur + tj^) + -y [£ r (n>) f + ri^rv)^ 

^ w 

(H y ) 4 = + VS,) - M — 


(H^)5 = 0 


(B. 13a) 
(B.13b) 


(B.13c). 


(B.13d) 

(B.13e) 
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Performing the linearization, the Jacobian coefficient matrix dEvJdQ is 


where 



a xx =(2 P z + X^ 2 + ^ 2 

a„ = rt 2 + (2 M + X)Z r 2 

*zz = x + M^r 2 ' 
a xr = (M 

=7-^r 


«0 : 


J_ t 2 
r jr 

k (Sx 2 + tf) 





(B.14) 
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The Jacobian coefficient matrix for the remaining non-cross derivative viscous terms, dF vJdQ, has the 

A A 

same form as dFyJdQ, but with £ replaced by rj. 

A A A 

And finally, linearizing H Vf the Jacobian coefficient matrix dHvfdQ is 


where 





+ [ 2 m + A(i,r e + rtfj] TT + 



Proteus 2-D Analysis Description 


B. Axisymmetric Analysis 69 



[2 



/ W \ 


VP/ 

/ 1 \ P 1 . d ( 

' 1 \ 

\T)--T + m '-d^\ 

, P ) 


B.2.3 Equation Of State 

The equation of state given in Section 43 must be modified slightly to add the swirl velocity w. Thus, 

p={y— l)j^£ r -yp(w 2 + v 2 + w 2 )J (B.16) 

or, in terms of temperature, 

r=i[-^-|(^ + v 2 + w 2 )] (bit) 

The derivatives arising from the linearization are the same as those presented in Section 4.3, except for 

-?£- = L^-(u 2 +v 2 + w 2 ) (B.18a) 


dp y - 1 / 2 , 2 , 2\ 


dp 2 

dp 


= ~{y - i) w 


-£--4 


dp Cy [_ P V 

dT w 

d(pw) c yP 


If constant stagnation enthalpy can be assumed, the appropriate equation of state is 

p = ^—y — (“ 2 + y2 + ™ 2 )J 


and the temperature becomes 


T= ih-T (u2+v2+w2) l 


A gain , the derivatives arising from the linearization are the same as in Section 4.3, except for 

f- i f L k+|(« 2 +' 2 +- 2 )] 


dp y - 1 


(B.21a) 

(B.21b) 
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37 


w 


d(pw) c pP 

B.2.4 Linearized Governing Equation 

The linearized form of equation (B.8) can now be written as 


a 0, At 1 
AQ + 1 1 


(B.21c) 

(B.21d) 



l+0o r 


_0iAt_ j_ ) _a 
1+02 r \ 


At 1 ( d ( rF ) d ( r h At 1 f d ( rF V i ) a N 

■ + — 5 +H| + , , „ — | Ti r + n V 


AQ > = 


1 + 0, M d( dr. 


l+0o r 


d( 


dr, 


(1+0 3 )At [ ( d(rE y J d(rFy 2 ) \ 0 3 At ) ( d(rE Vi ) d(rFy 2 ) 

+ T 7 . ~ I TT + - I ■ I ~7 *■ I . + 


n - 1 


l+0 2 r 
0 2 




dr, 


l + 0o r 


dt 


dr, 


A Q" ~ 1 + O (fii - y - 0 2 )(At) 2 , (0 3 - 0,)(At) 2 , (At) 3 J 


(B.22) 


1 + 0 2 

B.3 SOLUTION PROCEDURE 

Letting LHS(B.22) represent the left hand side of equation (B.22), we can write 


0 ,At 1 
LHS(B.22) = K + F 


d I dE 

nr lr ^-- r 


dE v i \ a ( aF d?v \ 


dQ dQ 




dQ dQ 


AQ n 


(B.23a) 


where 


K = l + 


1 / dH 


\+e 2 r 


(B.23b) 


dQ dQ J 

and I represents the identity matrix. The term in braces in equation (B.23a) can be factored to give 


LHS(B.22) = 


K+ I^±j r ( r jL_ r dEv ' 


\+0 2 r at. 


( e _Aj \ 2 i ir 1 

V 1+& 2 ) r 2 


dQ dQ 


, + K ->^_J -J-( r jL_ r dTv ' 


dE v 


' dQ dQ I dr< 


1 + 0 2 r dr, y aQ 
' * dEy 


AQ 


dQ dQ 


AQ 


(B.24) 


The last term represents the splitting error. 
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Equation (B.22) can thus be rewritten in spatially factored form, and, neglecting the temporal truncation 
and splitting error terms, becomes 


K + - 



,«L_, 

*\ 

dQ 

V 

Sir E) 

d(r F) 

a? 

dr} 

/ ^E K2 ) 

' l 35 + 


d£. 


dQ 


(K 'O' 


K + - 


+ A" a, l( **£+*!£+& 

1 + 0 2 r ^ dl drj J + 1 + $ 2 r ^ dl drj 



d¥ 

r Sr, \ 

V * 

A \ 

^f K j ) 

+ h k 

drj 


d¥ v 


dQ 


AQ = 


drf 


1 r 


dl 


drj 


+ TTe;^ n ~ l 


(B.25) 


Using the procedure of Douglas and Gunn (1964), as written by Briley and McDonald (1977), equation 
(B.25) can be split into the following two-sweep sequence. 

Sweep 1 (£ direction) 


K+ ^‘ AT 

i a 

(,jL 

_ r <\ 

K+ l+0 2 

r a? 1 

( dQ 

dQ J 


AQ = 


d(rE v ) d(r¥ v ) 
d \~ + drj 


Ax ! ( d(r E) , d(r¥) i , Ax ±( 

1 + e 2 r y dl + drj + + 1 + e 2 r y 

(1+<9 3 )Ax j ( d(rt V2 ) d(r¥ v }\ n B . jAt , / a(rF^) \" * ^ 

+ 1 + 6 2 r I 3{ + dr, J 1 + 0, *■ l 3{ + dr, J + \ + 0 2 


AQ 


(B.26a) 


Sweep 2 (r? direction) 


K+ 4^_J.J_[ r JL- r 


1+^2 r dtj 


dQ 


d¥ v 


dQ 


Or, expanding K and rearranging, 
Sweep 1 (£ direction) 


AQ" = K"AQ* 


(B.26b) 


1 d ( dZ \ 

'+*2 r * [ VaQ )' 

Ax 1 / 

r V 


] d 

1 + 0 2 r dl 


aE v 


AQ 


d(r E) d(r¥) 

-ic • a__ ' 


A )'t^ 


a(/-E K[ ) 5(r F y t ) 


+ 


*1** ! / d H dH v 


1 + 0 2 r 


dQ dQ 


d( 


drj 


+ H* 


(1+0 3 )At 1 f d(r£y 2 ) d(r F V J \ 0 3 Ax , / d(r F^) 

"t" i-l -it “F I -i./irl ^ t ~F 


<3? 




l+0 2 ' 


dl 


dr} 


1 + 0 2 


2 AQ"" 1 


AQ = 


(B.27a) 
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Sweep 2 (n direction) 


a 0iA? 1 A 
AQ + - 1 1 0 


l + 0 2 ^ dr} 


r I -#• | AQ" 
5Q 


_£iA^_i__a_ 
1 + 0 2 r &V 


df v 


dQ 


AQ" 


1+0, r 


, *i At 1 / dH 5 H k \ An 
+ TTT- TI — 7 7- AQ = 


dQ dQ 


a . 0 ,At 1 / aH 0Hv \ a , 

AQ +TTTT "a“ — UQ 


(B.27b) 


1 4* 0, ^ 


dQ dQ 


Applying the spatial differencing formulas of Section 5.0 results in 
Sweep 1 (£ direction) 


AQ, + 


0i At 


(1 + 0 2 )2A£ 

0jAt ! 


dE 

dQ 


AQ, + 1“ [r 


3E 

dQ 


AQ t -i 


(1 + 0 2 )2(A£y 

ft At 1 ( dH d U v 
dQ dQ 


l +0, ' 


T [ft - ifi - 1 + '.-/>V- AQ<- 1 - ft- 1 /- 1 + ft/ + i/ + i)VAQ; + ft/ + r 4 . + 1 / + ,)V + i A Q*+ 1 ] 
A Q' = - TTV T IV => + V F) + A} 1 + T IV W + W Fk,) + H„]" 


+ ILt^l i. [ 5{ < r E, 2 ) + <5,(r -y&4|V E^) + <5,(r F^]"' 1 +T ^- AQ"' 1 


(B.28a) 


Sweep 2 (>7 direction) 


.\n n dlAr 1 
” (1 + 0 2 )2Af? ' 


aF 

dQ 


(1 + 0 2 )2(At/)‘ 

ftA? 1 ( an 


1 + 0 2 ' 


1 a4; -'l 

y+i x / y-i J 

f [_ {r j r jfj)X- i A ft 1 — 1 + T j+ifj+ \ )V A 9/ + M + r j+ifj + 1 >V+ i A ft 1 ] 

V - 

I AQ (B.28b) 


3Q 3Q 


,A" „A‘ . e i Ax 1 / 5H 

AQ =AQ + . — — — 

1 + \ 3Q 3Q 


These equations are solved using the same matrix inversion procedure described in Section 7.2. 
B.4 CHIEN k-e TURBULENCE MODEL 

The axisymmetric k-e equations for the Chien model can be written using vector notation as 

d(rW) d(r¥) d(rG) 


dt 




dx 




dr 


= r{S + T) 


(B.29) 


where W, F, G, S, and T are the same as the corresponding terms in the two-dimensional planar equations 
with the coordinate y replaced by r , and 


-((.t) ; + (t) ! + (T) 2 ]-|(f + t + r) 2 

+ (£ + £) !+ (£) 2+ Km] 2 


(B.30a) 
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(B.30b) 


?2 = 


du dv_ _v 

dx dr r 


The analysis for the axisymmetric k-t equations is the same as for the two-dimensional planar k-t equations, 
described in the main body of this report. 
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ity model. The thin-layer or Euler equations may also be solved. The energy equation may be eliminated by the 
assumption of constant total enthalpy. Explicit and implicit artificial viscosity may be used. Several time step options 
are available for convergence acceleration. The documentation is divided into three volumes. This is the Analysis 
Description, and presents the equations and solution procedure. It describes in detail the governing equations, the 
turbulence model, the linearization of the equations and boundary conditions, the time and space differencing formu- 
las, the ADI solution procedure, and the artificial viscosity models. 
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